diff --git a/Mesh/highOrderSmoother.cpp b/Mesh/highOrderSmoother.cpp
index 5b3d6fa214279f726ad0273191749275e84f0a29..85c1e8850cd88ab2b6ea96ec6a4c648a22181498 100644
--- a/Mesh/highOrderSmoother.cpp
+++ b/Mesh/highOrderSmoother.cpp
@@ -100,9 +100,9 @@ void highOrderSmoother::updateTargetLocation(MVertex*v, const SPoint3 &p3,
struct p2data{
GFace *gf;
- MTriangle *t1,*t2;
+ MTriangle *t1, *t2;
MVertex *n12;
- fullMatrix<double> *m1,*m2;
+ fullMatrix<double> *m1, *m2;
highOrderSmoother *s;
p2data(GFace *_gf, MTriangle *_t1, MTriangle *_t2, MVertex *_n12,
highOrderSmoother *_s)
@@ -111,16 +111,18 @@ struct p2data{
elasticityTerm el(0, 1.e3, .3333,1);
s->moveToStraightSidedLocation(t1);
s->moveToStraightSidedLocation(t2);
- m1 = new fullMatrix<double>(3 * t1->getNumVertices(),
- 3 * t1->getNumVertices());
- m2 = new fullMatrix<double>(3 * t2->getNumVertices(),
- 3 * t2->getNumVertices());
- el.elementMatrix(t1,*m1);
- el.elementMatrix(t2,*m2);
+ m1 = new fullMatrix<double>(3 * t1->getNumVertices(),
+ 3 * t1->getNumVertices());
+ m2 = new fullMatrix<double>(3 * t2->getNumVertices(),
+ 3 * t2->getNumVertices());
+ SElement se1(t1), se2(t2);
+ el.elementMatrix(&se1, *m1);
+ el.elementMatrix(&se2, *m2);
s->moveToTargetLocation(t1);
s->moveToTargetLocation(t2);
}
- ~p2data(){
+ ~p2data()
+ {
delete m1;
delete m2;
}
@@ -143,12 +145,14 @@ struct pNdata{
3 * t1->getNumVertices());
m2 = new fullMatrix<double>(3 * t2->getNumVertices(),
3 * t2->getNumVertices());
- el.elementMatrix(t1,*m1);
- el.elementMatrix(t2,*m2);
+ SElement se1(t1), se2(t2);
+ el.elementMatrix(&se1, *m1);
+ el.elementMatrix(&se2, *m2);
s->moveToTargetLocation(t1);
s->moveToTargetLocation(t2);
}
- ~pNdata(){
+ ~pNdata()
+ {
delete m1;
delete m2;
}
@@ -551,7 +555,7 @@ double highOrderSmoother::smooth_metric_(std::vector<MElement*> & v,
if (myAssembler.sizeOfR()){
for (unsigned int i = 0; i < v.size(); i++){
- MElement *e = v[i];
+ MElement *e = v[i];
int nbNodes = e->getNumVertices();
const int n2 = 2 * nbNodes;
const int n3 = 3 * nbNodes;
@@ -564,21 +568,18 @@ double highOrderSmoother::smooth_metric_(std::vector<MElement*> & v,
fullVector<double> R2(n2);
fullMatrix<double> J23K33(n2, n3);
K33.set_all(0.0);
- El.elementMatrix(e, K33);
+ SElement se(e);
+ El.elementMatrix(&se, K33);
computeMetricVector(gf, e, J32, J23, D3);
J23K33.gemm(J23, K33, 1, 0);
K22.gemm(J23K33, J32, 1, 0);
J23K33.mult(D3, R2);
for (int j = 0; j < n2; j++){
- MVertex *vR;
- int iCompR, iFieldR;
- Dof RDOF = El.getLocalDofR(e, j);
- myAssembler.assemble(RDOF,-R2(j));
+ Dof RDOF = El.getLocalDofR(&se, j);
+ myAssembler.assemble(RDOF, -R2(j));
for (int k = 0; k < n2; k++){
- MVertex *vC;
- int iCompC, iFieldC;
- Dof CDOF = El.getLocalDofC(e, k);
- myAssembler.assemble(RDOF,CDOF,K22(j, k));
+ Dof CDOF = El.getLocalDofC(&se, k);
+ myAssembler.assemble(RDOF, CDOF, K22(j, k));
}
}
}
@@ -697,9 +698,10 @@ void highOrderSmoother::smooth(std::vector<MElement*> &all)
// assembly of the elasticity term on the
// set of elements
- for (int i=0;i<v.size();i++)
- El.addToMatrix(myAssembler, v[i]);
-
+ for (int i = 0; i < v.size(); i++){
+ SElement se(v[i]);
+ El.addToMatrix(myAssembler, &se);
+ }
// solve the system
lsys->systemSolve();
}
diff --git a/Solver/SElement.h b/Solver/SElement.h
new file mode 100644
index 0000000000000000000000000000000000000000..b9a658f4bc7ff6c6e2a54b3445b4efdc08f9db46
--- /dev/null
+++ b/Solver/SElement.h
@@ -0,0 +1,26 @@
+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef _SELEMENT_H_
+#define _SELEMENT_H_
+
+#include "MElement.h"
+
+// A solver element.
+
+class SElement
+{
+ private:
+ // the underlying mesh element
+ MElement *_e;
+ // store discrete function space and other data here
+ // ...
+ public:
+ SElement(MElement *e) : _e(e) {}
+ ~SElement(){}
+ MElement *getMeshElement() const { return _e; }
+};
+
+#endif
diff --git a/Solver/elasticityTerm.cpp b/Solver/elasticityTerm.cpp
index 1d10aea44c3a704955f371752bb7caacdcae7d09..d67462bafd2d1ee25632edad7bc6cc43c4982f23 100644
--- a/Solver/elasticityTerm.cpp
+++ b/Solver/elasticityTerm.cpp
@@ -6,8 +6,9 @@
#include "elasticityTerm.h"
#include "Numeric.h"
-void elasticityTerm::elementMatrix(MElement *e, fullMatrix<double> &m) const
+void elasticityTerm::elementMatrix(SElement *se, fullMatrix<double> &m) const
{
+ MElement *e = se->getMeshElement();
int nbNodes = e->getNumVertices();
int integrationOrder = 3 * (e->getPolynomialOrder() - 1) ;
int npts;
@@ -75,8 +76,9 @@ void elasticityTerm::elementMatrix(MElement *e, fullMatrix<double> &m) const
}
}
-void elasticityTerm::elementVector(MElement *e, fullVector<double> &m) const
+void elasticityTerm::elementVector(SElement *se, fullVector<double> &m) const
{
+ MElement *e = se->getMeshElement();
int nbNodes = e->getNumVertices();
int integrationOrder = 2 * e->getPolynomialOrder();
int npts;
diff --git a/Solver/elasticityTerm.h b/Solver/elasticityTerm.h
index a08b0cbf463727bf8d1359e80260aead91dce494..2f3dc92f3fa8c5ae061e832cf5bce6748b022f8d 100644
--- a/Solver/elasticityTerm.h
+++ b/Solver/elasticityTerm.h
@@ -9,7 +9,7 @@
#include "femTerm.h"
#include "Gmsh.h"
#include "GModel.h"
-#include "MElement.h"
+#include "SElement.h"
#include "fullMatrix.h"
class elasticityTerm : public femTerm<double, double> {
@@ -19,12 +19,19 @@ class elasticityTerm : public femTerm<double, double> {
SVector3 _volumeForce;
public:
// element matrix size : 3 dofs per vertex
- virtual int sizeOfR(MElement *e) const { return 3 * e->getNumVertices(); }
- virtual int sizeOfC(MElement *e) const { return 3 * e->getNumVertices(); }
+ virtual int sizeOfR(SElement *se) const
+ {
+ return 3 * se->getMeshElement()->getNumVertices();
+ }
+ virtual int sizeOfC(SElement *se) const
+ {
+ return 3 * se->getMeshElement()->getNumVertices();
+ }
// order dofs in the local matrix :
// dx1, dx2 ... dxn, dy1, ..., dyn, dz1, ... dzn
- Dof getLocalDofR(MElement *e, int iRow) const
+ Dof getLocalDofR(SElement *se, int iRow) const
{
+ MElement *e = se->getMeshElement();
int iCompR = iRow / e->getNumVertices();
int ithLocalVertex = iRow % e->getNumVertices();
return Dof(e->getVertex(ithLocalVertex)->getNum(),
@@ -34,8 +41,8 @@ class elasticityTerm : public femTerm<double, double> {
elasticityTerm(GModel *gm, double E, double nu, int iField) :
femTerm<double, double>(gm), _E(E), _nu(nu), _iField(iField) {}
void setVector(const SVector3 &f) { _volumeForce = f; }
- void elementMatrix(MElement *e, fullMatrix<double> &m) const;
- void elementVector(MElement *e, fullVector<double> &m) const;
+ void elementMatrix(SElement *se, fullMatrix<double> &m) const;
+ void elementVector(SElement *se, fullVector<double> &m) const;
};
#endif
diff --git a/Solver/femTerm.h b/Solver/femTerm.h
index b221f886fadfe94cf2720ce6f5c9cf45a5b20d0d..29372dff928a9fca050181b7784b5d2bdb8a7abb 100644
--- a/Solver/femTerm.h
+++ b/Solver/femTerm.h
@@ -13,7 +13,7 @@
#include "simpleFunction.h"
#include "dofManager.h"
#include "GModel.h"
-#include "MElement.h"
+#include "SElement.h"
// a nodal finite element term : variables are always defined at nodes
// of the mesh
@@ -25,20 +25,20 @@ class femTerm {
femTerm(GModel *gm) : _gm(gm) {}
virtual ~femTerm(){}
// return the number of columns of the element matrix
- virtual int sizeOfC(MElement*) const = 0;
+ virtual int sizeOfC(SElement *se) const = 0;
// return the number of rows of the element matrix
- virtual int sizeOfR(MElement*) const = 0;
+ virtual int sizeOfR(SElement *se) const = 0;
// in a given element, return the dof associated to a given row (column)
// of the local element matrix
- virtual Dof getLocalDofR(MElement *e, int iRow) const = 0;
+ virtual Dof getLocalDofR(SElement *se, int iRow) const = 0;
// default behavior: symmetric
- virtual Dof getLocalDofC(MElement *e, int iCol) const
+ virtual Dof getLocalDofC(SElement *se, int iCol) const
{
- return getLocalDofR(e, iCol);
+ return getLocalDofR(se, iCol);
}
// compute the elementary matrix
- virtual void elementMatrix(MElement *e, fullMatrix<dataMat> &m) const = 0;
- virtual void elementVector(MElement *e, fullVector<dataVec> &m) const
+ virtual void elementMatrix(SElement *se, fullMatrix<dataMat> &m) const = 0;
+ virtual void elementVector(SElement *se, fullVector<dataVec> &m) const
{
m.scale(0.0);
}
@@ -47,29 +47,29 @@ class femTerm {
void addToMatrix(dofManager<dataVec, dataMat> &dm, GEntity *ge) const
{
for(unsigned int i = 0; i < ge->getNumMeshElements(); i++){
- MElement *e = ge->getMeshElement(i);
- addToMatrix(dm, e);
+ SElement se(ge->getMeshElement(i));
+ addToMatrix(dm, &se);
}
}
// add the contribution from a single element to the dof manager
- void addToMatrix(dofManager<dataVec, dataMat> &dm, MElement *e) const
+ void addToMatrix(dofManager<dataVec, dataMat> &dm, SElement *se) const
{
- const int nbR = sizeOfR(e);
- const int nbC = sizeOfC(e);
+ const int nbR = sizeOfR(se);
+ const int nbC = sizeOfC(se);
fullMatrix<dataMat> localMatrix(nbR, nbC);
- elementMatrix(e, localMatrix);
- addToMatrix(dm, localMatrix, e);
+ elementMatrix(se, localMatrix);
+ addToMatrix(dm, localMatrix, se);
}
void addToMatrix(dofManager<dataVec, dataMat> &dm,
fullMatrix<dataMat> &localMatrix,
- MElement *e) const
+ SElement *se) const
{
const int nbR = localMatrix.size1();
const int nbC = localMatrix.size2();
for (int j = 0; j < nbR; j++){
- Dof R = getLocalDofR(e, j);
+ Dof R = getLocalDofR(se, j);
for (int k = 0; k < nbC; k++){
- Dof C = getLocalDofC(e, k);
+ Dof C = getLocalDofC(se, k);
dm.assemble(R, C, localMatrix(j, k));
}
}
@@ -120,15 +120,15 @@ class femTerm {
}
}
}
- void addToRightHandSide(dofManager<dataVec,dataMat> &dm, GEntity *ge) const
+ void addToRightHandSide(dofManager<dataVec, dataMat> &dm, GEntity *ge) const
{
for(unsigned int i = 0; i < ge->getNumMeshElements(); i++){
- MElement *e = ge->getMeshElement(i);
- int nbR = sizeOfR(e);
+ SElement se(ge->getMeshElement(i));
+ int nbR = sizeOfR(&se);
fullVector<dataVec> V(nbR);
- elementVector(e, V);
+ elementVector(&se, V);
// assembly
- for (int j = 0; j < nbR; j++) dm.assemble(getLocalDofR(e, j), V(j));
+ for (int j = 0; j < nbR; j++) dm.assemble(getLocalDofR(&se, j), V(j));
}
}
};
diff --git a/Solver/helmholtzTerm.h b/Solver/helmholtzTerm.h
index 01b82a60f13214c6e5c22ed151b20771aed4e001..3cf3e9e4f3f19fb7db0c784bfb0ca50e1dc9b96e 100644
--- a/Solver/helmholtzTerm.h
+++ b/Solver/helmholtzTerm.h
@@ -11,7 +11,7 @@
#include "simpleFunction.h"
#include "Gmsh.h"
#include "GModel.h"
-#include "MElement.h"
+#include "SElement.h"
#include "fullMatrix.h"
#include "Numeric.h"
@@ -28,18 +28,25 @@ class helmoltzTerm : public femTerm<scalar, scalar> {
: femTerm<scalar,scalar>(gm), _k(k), _a(a), _iFieldR(iFieldR),
_iFieldC(iFieldC) {}
// one dof per vertex (nodal fem)
- virtual int sizeOfR(MElement *e) const { return e->getNumVertices(); }
- virtual int sizeOfC(MElement *e) const { return e->getNumVertices(); }
- Dof getLocalDofR(MElement *e, int iRow) const
+ virtual int sizeOfR(SElement *se) const
+ {
+ return se->getMeshElement()->getNumVertices();
+ }
+ virtual int sizeOfC(SElement *se) const
+ {
+ return se->getMeshElement()->getNumVertices();
+ }
+ Dof getLocalDofR(SElement *se, int iRow) const
{
- return Dof(e->getVertex(iRow)->getNum(), _iFieldR);
+ return Dof(se->getMeshElement()->getVertex(iRow)->getNum(), _iFieldR);
}
- Dof getLocalDofC(MElement *e, int iRow) const
+ Dof getLocalDofC(SElement *se, int iRow) const
{
- return Dof(e->getVertex(iRow)->getNum(), _iFieldC);
+ return Dof(se->getMeshElement()->getVertex(iRow)->getNum(), _iFieldC);
}
- virtual void elementMatrix(MElement *e, fullMatrix<scalar> &m) const
+ virtual void elementMatrix(SElement *se, fullMatrix<scalar> &m) const
{
+ MElement *e = se->getMeshElement();
// compute integration rule
const int integrationOrder = (_a) ? 2 * e->getPolynomialOrder() :
2 * (e->getPolynomialOrder() - 1);