From 19c91ecff62d8399a1ec4a0dd08e5f6535faf084 Mon Sep 17 00:00:00 2001
From: Gaetan Bricteux <gaetan.bricteux@uclouvain.be>
Date: Wed, 27 Oct 2010 15:39:16 +0000
Subject: [PATCH] fix integration of polys + pp
---
Geo/MElement.cpp | 28 ++--
Geo/MElement.h | 3 +
Geo/MElementCut.cpp | 176 ++++++++-------------
Geo/MElementCut.h | 74 +++++----
Mesh/HighOrder.cpp | 8 +-
Solver/FuncGradDisc.h | 68 ++++----
Solver/SElement.cpp | 4 +-
Solver/convexCombinationTerm.h | 24 +--
Solver/crossConfTerm.h | 33 ++--
Solver/diagBCTerm.h | 19 +--
Solver/distanceTerm.h | 3 +-
Solver/elasticityTerm.cpp | 44 +++---
Solver/elasticityTerm.h | 16 +-
Solver/functionSpace.cpp | 6 +-
Solver/functionSpace.h | 237 +++++++++++++---------------
Solver/helmholtzTerm.h | 18 +--
Solver/laplaceTerm.h | 18 +--
Solver/orthogonalTerm.h | 41 +++--
Solver/terms.h | 273 +++++++++++++++------------------
19 files changed, 502 insertions(+), 591 deletions(-)
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 1166d924f0..d8416bafe5 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -244,8 +244,8 @@ double MElement::getJacobian(double u, double v, double w, double jac[3][3])
double gsf[256][3];
getGradShapeFunctions(u, v, w, gsf);
- for (int i = 0; i < getNumVertices(); i++) {
- const MVertex* v = getVertex(i);
+ for (int i = 0; i < getNumShapeFunctions(); i++) {
+ const MVertex *v = getShapeFunctionNode(i);
double* gg = gsf[i];
for (int j = 0; j < 3; j++) {
jac[j][0] += v->x() * gg[j];
@@ -269,8 +269,8 @@ double MElement::getJacobian(const fullMatrix<double> &gsf, double jac[3][3])
jac[1][0] = jac[1][1] = jac[1][2] = 0.;
jac[2][0] = jac[2][1] = jac[2][2] = 0.;
- for (int i = 0; i < getNumVertices(); i++) {
- const MVertex* v = getVertex(i);
+ for (int i = 0; i < getNumShapeFunctions(); i++) {
+ const MVertex *v = getShapeFunctionNode(i);
for (int j = 0; j < 3; j++) {
jac[j][0] += v->x() * gsf(i, j);
jac[j][1] += v->y() * gsf(i, j);
@@ -289,8 +289,8 @@ double MElement::getPrimaryJacobian(double u, double v, double w, double jac[3][
double gsf[256][3];
getGradShapeFunctions(u, v, w, gsf, 1);
- for(int i = 0; i < getNumPrimaryVertices(); i++) {
- const MVertex* v = getVertex(i);
+ for(int i = 0; i < getNumPrimaryShapeFunctions(); i++) {
+ const MVertex *v = getShapeFunctionNode(i);
double* gg = gsf[i];
for (int j = 0; j < 3; j++) {
jac[j][0] += v->x() * gg[j];
@@ -307,8 +307,8 @@ void MElement::pnt(double u, double v, double w, SPoint3 &p)
double x = 0., y = 0., z = 0.;
double sf[256];
getShapeFunctions(u, v, w, sf);
- for (int j = 0; j < getNumVertices(); j++) {
- const MVertex* v = getVertex(j);
+ for (int j = 0; j < getNumShapeFunctions(); j++) {
+ const MVertex *v = getShapeFunctionNode(j);
x += sf[j] * v->x();
y += sf[j] * v->y();
z += sf[j] * v->z();
@@ -321,8 +321,8 @@ void MElement::primaryPnt(double u, double v, double w, SPoint3 &p)
double x = 0., y = 0., z = 0.;
double sf[256];
getShapeFunctions(u, v, w, sf, 1);
- for (int j = 0; j < getNumPrimaryVertices(); j++) {
- const MVertex* v = getVertex(j);
+ for (int j = 0; j < getNumPrimaryShapeFunctions(); j++) {
+ const MVertex *v = getShapeFunctionNode(j);
x += sf[j] * v->x();
y += sf[j] * v->y();
z += sf[j] * v->z();
@@ -345,8 +345,8 @@ void MElement::xyz2uvw(double xyz[3], double uvw[3])
double xn = 0., yn = 0., zn = 0.;
double sf[256];
getShapeFunctions(uvw[0], uvw[1], uvw[2], sf);
- for (int i = 0; i < getNumVertices(); i++) {
- MVertex *v = getVertex(i);
+ for (int i = 0; i < getNumShapeFunctions(); i++) {
+ MVertex *v = getShapeFunctionNode(i);
xn += v->x() * sf[i];
yn += v->y() * sf[i];
zn += v->z() * sf[i];
@@ -396,7 +396,7 @@ double MElement::interpolate(double val[], double u, double v, double w, int str
int j = 0;
double sf[256];
getShapeFunctions(u, v, w, sf, order);
- for(int i = 0; i < getNumVertices(); i++){
+ for(int i = 0; i < getNumShapeFunctions(); i++){
sum += val[j] * sf[i];
j += stride;
}
@@ -410,7 +410,7 @@ void MElement::interpolateGrad(double val[], double u, double v, double w, doubl
int j = 0;
double gsf[256][3];
getGradShapeFunctions(u, v, w, gsf, order);
- for(int i = 0; i < getNumVertices(); i++){
+ for(int i = 0; i < getNumShapeFunctions(); i++){
dfdu[0] += val[j] * gsf[i][0];
dfdu[1] += val[j] * gsf[i][1];
dfdu[2] += val[j] * gsf[i][2];
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 47e4c16e8b..f758c8c166 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -233,6 +233,9 @@ class MElement
double getJacobian(double u, double v, double w, double jac[3][3]);
double getPrimaryJacobian(double u, double v, double w, double jac[3][3]);
double getJacobianDeterminant(double u, double v, double w);
+ virtual int getNumShapeFunctions(){ return getNumVertices(); }
+ virtual int getNumPrimaryShapeFunctions(){ return getNumPrimaryVertices(); }
+ virtual MVertex *getShapeFunctionNode(int i){ return getVertex(i); }
// get the point in cartesian coordinates corresponding to the point
// (u,v,w) in parametric coordinates
diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index 76b030f4ea..2725ef0588 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -77,7 +77,8 @@ void MPolyhedron::_init()
// innerVertices
for(unsigned int i = 0; i < _parts.size(); i++) {
for(int j = 0; j < 4; j++) {
- if(std::find(_vertices.begin(), _vertices.end(), _parts[i]->getVertex(j)) == _vertices.end())
+ if(std::find(_vertices.begin(), _vertices.end(), _parts[i]->getVertex(j)) ==
+ _vertices.end())
_innerVertices.push_back(_parts[i]->getVertex(j));
}
}
@@ -136,13 +137,13 @@ void MPolyhedron::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
const double u = ptsi[ip].pt[0];
const double v = ptsi[ip].pt[1];
const double w = ptsi[ip].pt[2];
- const double weight = ptsi[ip].weight;
- const double detJ = tt.getJacobian(u, v, w, jac);
SPoint3 p; tt.pnt(u, v, w, p);
_intpt[*npts + ip].pt[0] = p.x();
_intpt[*npts + ip].pt[1] = p.y();
_intpt[*npts + ip].pt[2] = p.z();
- _intpt[*npts + ip].weight = detJ * weight;
+ double partJac = _parts[i]->getJacobian(p.x(), p.y(), p.z(), jac);
+ double Jac = getJacobian(p.x(), p.y(), p.z(), jac);
+ _intpt[*npts + ip].weight = ptsi[ip].weight * partJac / Jac;
}
*npts += nptsi;
}
@@ -204,21 +205,23 @@ void MPolygon::_initVertices()
// innerVertices
for(unsigned int i = 0; i < _parts.size(); i++) {
for(int j = 0; j < 3; j++) {
- if(std::find(_vertices.begin(), _vertices.end(), _parts[i]->getVertex(j)) == _vertices.end())
+ if(std::find(_vertices.begin(), _vertices.end(), _parts[i]->getVertex(j)) ==
+ _vertices.end())
_innerVertices.push_back(_parts[i]->getVertex(j));
}
}
}
+
bool MPolygon::isInside(double u, double v, double w)
{
- if(!_orig) return false;
+ if(!getParent()) return false;
double uvw[3] = {u, v, w};
for(unsigned int i = 0; i < _parts.size(); i++) {
double v_uvw[3][3];
for(int j = 0; j < 3; j++) {
MVertex *vij = _parts[i]->getVertex(j);
double v_xyz[3] = {vij->x(), vij->y(), vij->z()};
- _orig->xyz2uvw(v_xyz, v_uvw[j]);
+ getParent()->xyz2uvw(v_xyz, v_uvw[j]);
}
MVertex v0(v_uvw[0][0], v_uvw[0][1], v_uvw[0][2]);
MVertex v1(v_uvw[1][0], v_uvw[1][1], v_uvw[1][2]);
@@ -236,7 +239,7 @@ void MPolygon::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
{
*npts = 0;
if(_intpt) delete [] _intpt;
- if(!_orig) return;
+ if(!getParent()) return;
double jac[3][3];
_intpt = new IntPt[getNGQTPts(pOrder) * _parts.size()];
for(unsigned int i = 0; i < _parts.size(); i++) {
@@ -247,7 +250,7 @@ void MPolygon::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
for(int j = 0; j < 3; j++) {
double xyz[3] = {_parts[i]->getVertex(j)->x(), _parts[i]->getVertex(j)->y(),
_parts[i]->getVertex(j)->z()};
- _orig->xyz2uvw(xyz, uvw[j]);
+ getParent()->xyz2uvw(xyz, uvw[j]);
}
MVertex v0(uvw[0][0], uvw[0][1], uvw[0][2]);
MVertex v1(uvw[1][0], uvw[1][1], uvw[1][2]);
@@ -257,13 +260,13 @@ void MPolygon::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
const double u = ptsi[ip].pt[0];
const double v = ptsi[ip].pt[1];
const double w = ptsi[ip].pt[2];
- const double weight = ptsi[ip].weight;
- const double detJ = tt.getJacobian(u, v, w, jac);
SPoint3 p; tt.pnt(u, v, w, p);
_intpt[*npts + ip].pt[0] = p.x();
_intpt[*npts + ip].pt[1] = p.y();
_intpt[*npts + ip].pt[2] = p.z();
- _intpt[*npts + ip].weight = detJ * weight;
+ double partJac = _parts[i]->getJacobian(p.x(), p.y(), p.z(), jac);
+ double Jac = getJacobian(p.x(), p.y(), p.z(), jac);
+ _intpt[*npts + ip].weight = ptsi[ip].weight * partJac / Jac;
}
*npts += nptsi;
}
@@ -299,74 +302,55 @@ void MLineChild::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
if(!_orig) return;
double jac[3][3];
_intpt = new IntPt[getNGQLPts(pOrder)];
-
int nptsi;
IntPt *ptsi;
double v_uvw[2][3];
-
-// -------------before--------------------//
-
-// for(int i = 0; i < 2; i++) {
-// MVertex *vi = getVertex(i);
-// double v_xyz[3] = {vi->x(), vi->y(), vi->z()};
-// _orig->xyz2uvw(v_xyz, v_uvw[i]);
-// }
-
-// -----------mich mach---------------------//
-
-
- MVertex *vo = _orig->getVertex(0);
- MVertex *vf = _orig->getVertex(1);
-
- SPoint3 P(vo->x(), vo->y(), vo->z());
- SPoint3 Q(vf->x(), vf->y(), vf->z());
-
- SPoint3 R;
-
- R = P + Q;
- R/=2;
-
- vo = getVertex(0);
- vf = getVertex(1);
-
- SPoint3 A(vo->x(), vo->y(), vo->z());
- SPoint3 B(vf->x(), vf->y(), vf->z());
-
- double lengthDemi = R.distance(Q);
-
- if (P.distance(A) < Q.distance(A)) {v_uvw[0][0] = - R.distance(A)/lengthDemi;v_uvw[0][1]=0;v_uvw[0][2]=0;}
- else {v_uvw[0][0] = R.distance(A)/lengthDemi;v_uvw[0][1]=0;v_uvw[0][2]=0;}
-
- if (P.distance(B) < Q.distance(B)) {v_uvw[1][0] = - R.distance(B)/lengthDemi;v_uvw[1][1]=0;v_uvw[1][2]=0;}
- else {v_uvw[1][0] = R.distance(B)/lengthDemi;v_uvw[1][1]=0;v_uvw[1][2]=0;}
-
-// ---------------------------------------//
-
+ for(int i = 0; i < 2; i++) {
+ MVertex *vi = getVertex(i);
+ double v_xyz[3] = {vi->x(), vi->y(), vi->z()};
+ _orig->xyz2uvw(v_xyz, v_uvw[i]);
+ }
MVertex v0(v_uvw[0][0], v_uvw[0][1], v_uvw[0][2]);
MVertex v1(v_uvw[1][0], v_uvw[1][1], v_uvw[1][2]);
-
MLine l(&v0, &v1);
l.getIntegrationPoints(pOrder, &nptsi, &ptsi);
-
for(int ip = 0; ip < nptsi; ip++){
const double u = ptsi[ip].pt[0];
const double v = ptsi[ip].pt[1];
const double w = ptsi[ip].pt[2];
- const double weight = ptsi[ip].weight;
- const double detJ = l.getJacobian(u, v, w, jac);
SPoint3 p; l.pnt(u, v, w, p);
_intpt[*npts + ip].pt[0] = p.x();
_intpt[*npts + ip].pt[1] = p.y();
_intpt[*npts + ip].pt[2] = p.z();
- _intpt[*npts + ip].weight = detJ * weight;
+ _intpt[*npts + ip].weight = ptsi[ip].weight;
}
*npts = nptsi;
*pts = _intpt;
-
}
//----------------------------------- MTriangleBorder ------------------------------
+bool MTriangleBorder::isInside(double u, double v, double w)
+{
+ if(!getParent()) return false;
+ double uvw[3] = {u, v, w};
+ double v_uvw[3][3];
+ for(int i = 0; i < 3; i++) {
+ MVertex *vi = getVertex(i);
+ double v_xyz[3] = {vi->x(), vi->y(), vi->z()};
+ getParent()->xyz2uvw(v_xyz, v_uvw[i]);
+ }
+ MVertex v0(v_uvw[0][0], v_uvw[0][1], v_uvw[0][2]);
+ MVertex v1(v_uvw[1][0], v_uvw[1][1], v_uvw[1][2]);
+ MVertex v2(v_uvw[2][0], v_uvw[2][1], v_uvw[2][2]);
+ MTriangle t(&v0, &v1, &v2);
+ double ksi[3];
+ t.xyz2uvw(uvw, ksi);
+ if(t.isInside(ksi[0], ksi[1], ksi[2]))
+ return true;
+ return false;
+}
+
void MTriangleBorder::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
{
*npts = 0;
@@ -391,66 +375,39 @@ void MTriangleBorder::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
const double u = ptsi[ip].pt[0];
const double v = ptsi[ip].pt[1];
const double w = ptsi[ip].pt[2];
- const double weight = ptsi[ip].weight;
- //const double detJ =
tt.getJacobian(u, v, w, jac);
SPoint3 p; tt.pnt(u, v, w, p);
_intpt[ip].pt[0] = p.x();
_intpt[ip].pt[1] = p.y();
_intpt[ip].pt[2] = p.z();
- _intpt[ip].weight = weight;//detJ * weight;
+ _intpt[ip].weight = ptsi[ip].weight;
}
*npts = nptsi;
*pts = _intpt;
}
-//----------------------------------- MPolygonBorder ------------------------------
+//-------------------------------------- MLineBorder ------------------------------
-void MPolygonBorder::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
+bool MLineBorder::isInside(double u, double v, double w)
{
- *npts = 0;
- if(_intpt) delete [] _intpt;
- if(!getParent()) return;
- _intpt = new IntPt[getNGQTPts(pOrder) * _parts.size()];
- int nptsi;
- IntPt *ptsi;
-
- std::vector<MVertex*> verts;
- for(unsigned int i = 0; i < _parts.size(); i++) {
- double uvw[3][3];
- for(int j = 0; j < 3; j++) {
- MVertex *v = _parts[i]->getVertex(j);
- double xyz[3] = {v->x(), v->y(), v->z()};
- getParent()->xyz2uvw(xyz, uvw[j]);
- }
- verts.push_back(new MVertex(uvw[0][0], uvw[0][1], uvw[0][2]));
- verts.push_back(new MVertex(uvw[1][0], uvw[1][1], uvw[1][2]));
- verts.push_back(new MVertex(uvw[2][0], uvw[2][1], uvw[2][2]));
- }
- MPolygon pp(verts);
- pp.getIntegrationPoints(pOrder, &nptsi, &ptsi);
- double jac[3][3];
- for(int ip = 0; ip < nptsi; ip++){
- const double u = ptsi[ip].pt[0];
- const double v = ptsi[ip].pt[1];
- const double w = ptsi[ip].pt[2];
- const double weight = ptsi[ip].weight;
- //const double detJ =
- pp.getJacobian(u, v, w, jac);
- SPoint3 p; pp.pnt(u, v, w, p);
- _intpt[ip].pt[0] = p.x();
- _intpt[ip].pt[1] = p.y();
- _intpt[ip].pt[2] = p.z();
- _intpt[ip].weight = weight;//detJ * weight;
+ if(!getParent()) return false;
+ double uvw[3] = {u, v, w};
+ double v_uvw[2][3];
+ for(int i = 0; i < 2; i++) {
+ MVertex *vi = getVertex(i);
+ double v_xyz[3] = {vi->x(), vi->y(), vi->z()};
+ getParent()->xyz2uvw(v_xyz, v_uvw[i]);
}
- *npts = nptsi;
- *pts = _intpt;
- for(unsigned int i = 0; i < verts.size(); i++)
- delete verts[i];
+ MVertex v0(v_uvw[0][0], v_uvw[0][1], v_uvw[0][2]);
+ MVertex v1(v_uvw[1][0], v_uvw[1][1], v_uvw[1][2]);
+ MLine l(&v0, &v1);
+ double ksi[3];
+ l.xyz2uvw(uvw, ksi);
+ if(l.isInside(ksi[0], ksi[1], ksi[2]))
+ return true;
+ return false;
}
-//-------------------------------------- MLineBorder ------------------------------
-
void MLineBorder::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
{
*npts = 0;
@@ -474,14 +431,11 @@ void MLineBorder::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
const double u = ptsi[ip].pt[0];
const double v = ptsi[ip].pt[1];
const double w = ptsi[ip].pt[2];
- const double weight = ptsi[ip].weight;
- //const double detJ =
- ll.getJacobian(u, v, w, jac);
SPoint3 p; ll.pnt(u, v, w, p);
_intpt[ip].pt[0] = p.x();
_intpt[ip].pt[1] = p.y();
_intpt[ip].pt[2] = p.z();
- _intpt[ip].weight = weight;//detJ * weight;
+ _intpt[ip].weight = ptsi[ip].weight;
}
*npts = nptsi;
*pts = _intpt;
@@ -1235,7 +1189,8 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
cp.clear(); lines.clear(); triangles.clear(); quads.clear(); tetras.clear(); hexas.clear();
}
- /*for(int i = 0; i < 10; i++) {
+#if 0
+ for(int i = 0; i < 10; i++) {
printf(" - element type : %d\n", i);
for(std::map<int, std::vector<MElement*> >::iterator it = elements[i].begin();
it != elements[i].end(); it++){
@@ -1243,13 +1198,14 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
for(int j = 0; j < it->second.size(); j++){
MElement *e = it->second[j];
printf("element %d",e->getNum());
- if(e->getParent()) printf(" par=%d",e->getParent()->getNum());
+ if(e->getParent()) printf(" par=%d (%d)",e->getParent()->getNum(),e->ownsParent());
if(e->getDomain(0)) printf(" d0=%d",e->getDomain(0)->getNum());
if(e->getDomain(1)) printf(" d1=%d",e->getDomain(1)->getNum());
printf("\n");
}
}
- }printf("\n");*/
+ }printf("\n");
+#endif
for(newVerticesContainer::iterator it = newVertices.begin() ; it != newVertices.end(); ++it) {
vertexMap[(*it)->getNum()] = *it;
diff --git a/Geo/MElementCut.h b/Geo/MElementCut.h
index 137abb1615..e18bb4a8f7 100644
--- a/Geo/MElementCut.h
+++ b/Geo/MElementCut.h
@@ -116,13 +116,11 @@ class MPolyhedron : public MElement {
}
virtual const polynomialBasis* getFunctionSpace(int order=-1) const
{
- if(_orig) return _orig->getFunctionSpace(order);
- return 0;
+ return (_orig ? _orig->getFunctionSpace(order) : 0);
}
virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const
{
- if(_orig) return _orig->getJacobianFuncSpace(order);
- return 0;
+ return (_orig ? _orig->getJacobianFuncSpace(order) : 0);
}
virtual void getShapeFunctions(double u, double v, double w, double s[], int o)
{
@@ -132,12 +130,24 @@ class MPolyhedron : public MElement {
{
if(_orig) _orig->getGradShapeFunctions(u, v, w, s, o);
}
- // the parametric coordinates of the polyhedron are
- // the coordinates in the local parent element.
- virtual void xyz2uvw(double xyz[3], double uvw[3])
+ virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3], int o)
+ {
+ if(_orig) _orig->getHessShapeFunctions(u, v, w, s, o);
+ }
+ virtual int getNumShapeFunctions()
+ {
+ return (_orig ? _orig->getNumShapeFunctions() : 0);
+ }
+ virtual int getNumPrimaryShapeFunctions()
+ {
+ return (_orig ? _orig->getNumPrimaryShapeFunctions() : 0);
+ }
+ virtual MVertex *getShapeFunctionNode(int i)
{
- if(_orig) _orig->xyz2uvw(xyz,uvw);
+ return (_orig ? _orig->getShapeFunctionNode(i) : 0);
}
+ // the parametric coordinates of the polyhedron are
+ // the coordinates in the local parent element.
virtual bool isInside(double u, double v, double w);
virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
virtual MElement *getParent() const { return _orig; }
@@ -241,15 +251,18 @@ class MPolygon : public MElement {
_edges.clear();
_initVertices();
}
+ virtual MElement *getParent() const { return _orig; }
+ virtual void setParent(MElement *p, bool owner = false) { _orig = p; _owner = owner; }
+ virtual int getNumChildren() const { return _parts.size(); }
+ virtual MElement *getChild(int i) const { return _parts[i]; }
+ virtual bool ownsParent() const { return _owner; }
virtual const polynomialBasis* getFunctionSpace(int order=-1) const
{
- if(_orig) return _orig->getFunctionSpace(order);
- return 0;
+ return (_orig ? _orig->getFunctionSpace(order) : 0);
}
virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const
{
- if(_orig) return _orig->getJacobianFuncSpace(order);
- return 0;
+ return (_orig ? _orig->getJacobianFuncSpace(order) : 0);
}
virtual void getShapeFunctions(double u, double v, double w, double s[], int o)
{
@@ -259,19 +272,26 @@ class MPolygon : public MElement {
{
if(_orig) _orig->getGradShapeFunctions(u, v, w, s, o);
}
- // the parametric coordinates of the polygon are
- // the coordinates in the local parent element.
- virtual void xyz2uvw(double xyz[3], double uvw[3])
+ virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3], int o)
+ {
+ if(_orig) _orig->getHessShapeFunctions(u, v, w, s, o);
+ }
+ virtual int getNumShapeFunctions()
+ {
+ return (_orig ? _orig->getNumShapeFunctions() : 0);
+ }
+ virtual int getNumPrimaryShapeFunctions()
{
- if(_orig) _orig->xyz2uvw(xyz,uvw);
+ return (_orig ? _orig->getNumPrimaryShapeFunctions() : 0);
}
+ virtual MVertex *getShapeFunctionNode(int i)
+ {
+ return (_orig ? _orig->getShapeFunctionNode(i) : 0);
+ }
+ // the parametric coordinates of the polygon are
+ // the coordinates in the local parent element.
virtual bool isInside(double u, double v, double w);
virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
- virtual MElement *getParent() const { return _orig; }
- virtual void setParent(MElement *p, bool owner = false) { _orig = p; _owner = owner; }
- virtual int getNumChildren() const { return _parts.size(); }
- virtual MElement *getChild(int i) const { return _parts[i]; }
- virtual bool ownsParent() const { return _owner; }
virtual int getNumVerticesForMSH() {return _parts.size() * 3;}
virtual int *getVerticesIdForMSH()
{
@@ -320,12 +340,12 @@ class MLineChild : public MLine {
{
if(_orig) _orig->getGradShapeFunctions(u, v, w, s, o);
}
- // the parametric coordinates of the LineChildren are
- // the coordinates in the local parent element.
- virtual void xyz2uvw(double xyz[3], double uvw[3])
+ virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3], int o)
{
- if(_orig) _orig->xyz2uvw(xyz,uvw);
+ if(_orig) _orig->getHessShapeFunctions(u, v, w, s, o);
}
+ // the parametric coordinates of the LineChildren are
+ // the coordinates in the local parent element.
virtual bool isInside(double u, double v, double w);
virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
virtual MElement *getParent() const { return _orig; }
@@ -361,6 +381,7 @@ class MTriangleBorder : public MTriangle {
return NULL;
}
virtual int getTypeForMSH() const { return MSH_TRI_B; }
+ virtual bool isInside(double u, double v, double w);
// the integration points of the MTriangleBorder are in the parent element space
virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
};
@@ -391,8 +412,6 @@ class MPolygonBorder : public MPolygon {
return NULL;
}
virtual int getTypeForMSH() const { return MSH_POLYG_B; }
- // the integration points of the MPolygonBorder are in the parent element space
- virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
};
class MLineBorder : public MLine {
@@ -421,6 +440,7 @@ class MLineBorder : public MLine {
return NULL;
}
virtual int getTypeForMSH() const { return MSH_LIN_B; }
+ virtual bool isInside(double u, double v, double w);
// the integration points of the MLineBorder are in the parent element space
virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
};
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 3317c57fd0..6784b09e21 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -430,12 +430,10 @@ static void getFaceVertices(GFace *gf, MElement *incomplete, MElement *ele,
}
else{
double X(0), Y(0), Z(0), GUESS[2] = {0, 0};
- double sf[256];
- incomplete->getShapeFunctions(-1., -1., 0, sf);
- for (int j = 0; j < incomplete->getNumVertices(); j++)
+ double sf[256];
incomplete->getShapeFunctions(t1, t2, 0, sf);
- for (int j = 0; j < incomplete->getNumVertices(); j++){
- MVertex *vt = incomplete->getVertex(j);
+ for (int j = 0; j < incomplete->getNumShapeFunctions(); j++){
+ MVertex *vt = incomplete->getShapeFunctionNode(j);
X += sf[j] * vt->x();
Y += sf[j] * vt->y();
Z += sf[j] * vt->z();
diff --git a/Solver/FuncGradDisc.h b/Solver/FuncGradDisc.h
index 866fd50adc..a508c3aad2 100644
--- a/Solver/FuncGradDisc.h
+++ b/Solver/FuncGradDisc.h
@@ -36,24 +36,24 @@ class FuncGradDisc : public simpleFunctionOnElement<double> {
// --- F2 --- //
- MElement *e=getElement();
+ MElement *e = getElement();
SPoint3 p(x,y,z);
if (e->getParent()) e = e->getParent();
- int numV = e->getNumVertices();
double xyz[3] = {x,y,z};
double uvw[3];
e->xyz2uvw(xyz,uvw);
double val[30];
e->getShapeFunctions(uvw[0], uvw[1], uvw[2], val);
double f = 0;
- for (int i = 0;i<numV;i++)
+ for (int i = 0; i < e->getNumShapeFunctions(); i++)
{
+ MVertex *v = e->getShapeFunctionNode(i);
//std::cout<<"val[i] :" << val[i] << "\n";
- //std::cout<<"ls(i) :" << (*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z()) << "\n";
- f = f + std::abs((*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z()))*val[i];
+ //std::cout<<"ls(i) :" << (*_ls)(v->x(),v->y(),v->z()) << "\n";
+ f = f + std::abs((*_ls)(v->x(), v->y(), v->z())) * val[i];
}
- f = f - std::abs((*_ls)(x,y,z));
+ f = f - std::abs((*_ls)(x, y, z));
//std::cout<<"val f :" << f << "\n";
return f;
@@ -64,16 +64,16 @@ class FuncGradDisc : public simpleFunctionOnElement<double> {
// SPoint3 p(x,y,z);
// if (e->getParent()) e = e->getParent();
-// int numV = e->getNumVertices();
// double xyz[3] = {x,y,z};
// double uvw[3];
// e->xyz2uvw(xyz,uvw);
// double val[30];
// e->getShapeFunctions(uvw[0], uvw[1], uvw[2], val);
// double f = 0;
-// for (int i = 0;i<numV;i++)
+// for (int i = 0; i < e->getNumShapeFunctions(); i++)
// {
-// f = f + (*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z())*val[i];
+// MVertex *v = e-<getShapeFunctionNode(i);
+// f = f + (*_ls)(v->x(), v->y(), v->z()) * val[i];
// }
// f = std::abs(f);
// return f;
@@ -91,7 +91,6 @@ class FuncGradDisc : public simpleFunctionOnElement<double> {
MElement *e=getElement();
SPoint3 p(x,y,z);
if (e->getParent()) e = e->getParent();
- int numV = e->getNumVertices();
double xyz[3] = {x,y,z};
double uvw[3];
e->xyz2uvw(xyz,uvw);
@@ -112,33 +111,35 @@ class FuncGradDisc : public simpleFunctionOnElement<double> {
if ((*_ls)(x,y,z)>0)
{
- for (int i = 0;i<numV;i++)
+ for (int i = 0; i < e->getNumShapeFunctions(); i++)
{
dNdx = invjac[0][0] * gradsuvw[i][0] + invjac[0][1] * gradsuvw[i][1] + invjac[0][2] * gradsuvw[i][2];
dNdy = invjac[1][0] * gradsuvw[i][0] + invjac[1][1] * gradsuvw[i][1] + invjac[1][2] * gradsuvw[i][2];
dNdz = invjac[2][0] * gradsuvw[i][0] + invjac[2][1] * gradsuvw[i][1] + invjac[2][2] * gradsuvw[i][2];
- dfdx = dfdx + std::abs((*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z()))*dNdx ;
- dfdx = dfdx - (*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z())*dNdx;
- dfdy = dfdy + std::abs((*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z()))*dNdy ;
- dfdy = dfdy - (*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z())*dNdy;
- dfdz = dfdz + std::abs((*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z()))*dNdz ;
- dfdz = dfdz - (*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z())*dNdz;
+ MVertex *v = e->getShapeFunctionNode(i);
+ dfdx = dfdx + std::abs((*_ls)(v->x(), v->y(), v->z())) * dNdx ;
+ dfdx = dfdx - (*_ls)(v->x(), v->y(), v->z()) * dNdx;
+ dfdy = dfdy + std::abs((*_ls)(v->x(), v->y(), v->z())) * dNdy ;
+ dfdy = dfdy - (*_ls)(v->x(), v->y(), v->z()) * dNdy;
+ dfdz = dfdz + std::abs((*_ls)(v->x(), v->y(), v->z())) * dNdz ;
+ dfdz = dfdz - (*_ls)(v->x(), v->y(), v->z()) * dNdz;
}
}else
{
- for (int i = 0;i<numV;i++)
+ for (int i = 0; i < e->getNumShapeFunctions(); i++)
{
dNdx = invjac[0][0] * gradsuvw[i][0] + invjac[0][1] * gradsuvw[i][1] + invjac[0][2] * gradsuvw[i][2];
dNdy = invjac[1][0] * gradsuvw[i][0] + invjac[1][1] * gradsuvw[i][1] + invjac[1][2] * gradsuvw[i][2];
dNdz = invjac[2][0] * gradsuvw[i][0] + invjac[2][1] * gradsuvw[i][1] + invjac[2][2] * gradsuvw[i][2];
- dfdx = dfdx + std::abs((*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z()))*dNdx ;
- dfdx = dfdx + (*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z())*dNdx;
- dfdy = dfdy + std::abs((*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z()))*dNdy ;
- dfdy = dfdy + (*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z())*dNdy;
- dfdz = dfdz + std::abs((*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z()))*dNdz ;
- dfdz = dfdz + (*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z())*dNdz;
+ MVertex *v = e->getShapeFunctionNode(i);
+ dfdx = dfdx + std::abs((*_ls)(v->x(), v->y(), v->z())) * dNdx ;
+ dfdx = dfdx + (*_ls)(v->x(), v->y(), v->z()) * dNdx;
+ dfdy = dfdy + std::abs((*_ls)(v->x(), v->y(), v->z())) * dNdy ;
+ dfdy = dfdy + (*_ls)(v->x(), v->y(), v->z()) * dNdy;
+ dfdz = dfdz + std::abs((*_ls)(v->x(), v->y(), v->z())) * dNdz ;
+ dfdz = dfdz + (*_ls)(v->x(), v->y(), v->z()) * dNdz;
}
}
}
@@ -149,7 +150,6 @@ class FuncGradDisc : public simpleFunctionOnElement<double> {
//
// SPoint3 p(x,y,z);
// if (e->getParent()) e = e->getParent();
-// int numV = e->getNumVertices();
// double xyz[3] = {x,y,z};
// double uvw[3];
// e->xyz2uvw(xyz,uvw);
@@ -170,27 +170,29 @@ class FuncGradDisc : public simpleFunctionOnElement<double> {
//
// if ((*_ls)(x,y,z)>0)
// {
-// for (int i = 0;i<numV;i++)
+// for (int i = 0; i < e->getNumShapeFunctions(); i++)
// {
// dNdx = invjac[0][0] * gradsuvw[i][0] + invjac[0][1] * gradsuvw[i][1] + invjac[0][2] * gradsuvw[i][2];
// dNdy = invjac[1][0] * gradsuvw[i][0] + invjac[1][1] * gradsuvw[i][1] + invjac[1][2] * gradsuvw[i][2];
// dNdz = invjac[2][0] * gradsuvw[i][0] + invjac[2][1] * gradsuvw[i][1] + invjac[2][2] * gradsuvw[i][2];
//
-// dfdx = dfdx + (*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z())*dNdx;
-// dfdy = dfdy + (*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z())*dNdy;
-// dfdz = dfdz + (*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z())*dNdz;
+// MVertex *v = e->getShapeFunctionNode(i);
+// dfdx = dfdx + (*_ls)(v->x(), v->y(), v->z()) * dNdx;
+// dfdy = dfdy + (*_ls)(v->x(), v->y(), v->z()) * dNdy;
+// dfdz = dfdz + (*_ls)(v->x(), v->y(), v->z()) * dNdz;
// }
// }else
// {
-// for (int i = 0;i<numV;i++)
+// for (int i = 0; i < e->getNumShapeFunctions(); i++)
// {
// dNdx = invjac[0][0] * gradsuvw[i][0] + invjac[0][1] * gradsuvw[i][1] + invjac[0][2] * gradsuvw[i][2];
// dNdy = invjac[1][0] * gradsuvw[i][0] + invjac[1][1] * gradsuvw[i][1] + invjac[1][2] * gradsuvw[i][2];
// dNdz = invjac[2][0] * gradsuvw[i][0] + invjac[2][1] * gradsuvw[i][1] + invjac[2][2] * gradsuvw[i][2];
//
-// dfdx = dfdx - (*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z())*dNdx;
-// dfdy = dfdy - (*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z())*dNdy;
-// dfdz = dfdz - (*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z())*dNdz;
+// MVertex *v = e->getShapeFunctionNode(i);
+// dfdx = dfdx - (*_ls)(v->x(), v->y(), v->z()) * dNdx;
+// dfdy = dfdy - (*_ls)(v->x(), v->y(), v->z()) * dNdy;
+// dfdz = dfdz - (*_ls)(v->x(), v->y(), v->z()) * dNdz;
// }
// }
// }
diff --git a/Solver/SElement.cpp b/Solver/SElement.cpp
index 1c7a02354f..746c7d44c3 100644
--- a/Solver/SElement.cpp
+++ b/Solver/SElement.cpp
@@ -25,8 +25,8 @@ void SElement::gradNodalFunctions (double u, double v, double w, double invjac[3
double grads[256][3];
_e->getGradShapeFunctions(u, v, w, grads);
- int nbNodes = getNumNodalShapeFunctions();
- for (int j = 0; j < nbNodes; j++){
+ int nbSF = getNumNodalShapeFunctions();
+ for (int j = 0; j < nbSF; j++){
Grads[j][0] = invjac[0][0] * grads[j][0] + invjac[0][1] * grads[j][1] +
invjac[0][2] * grads[j][2];
Grads[j][1] = invjac[1][0] * grads[j][0] + invjac[1][1] * grads[j][1] +
diff --git a/Solver/convexCombinationTerm.h b/Solver/convexCombinationTerm.h
index e559fd3e21..b2e5b79e5c 100644
--- a/Solver/convexCombinationTerm.h
+++ b/Solver/convexCombinationTerm.h
@@ -25,15 +25,15 @@ class convexCombinationTerm : public femTerm<double> {
: femTerm<double>(gm), _k(k), _iField(iField) {}
virtual int sizeOfR(SElement *se) const
{
- return se->getMeshElement()->getNumVertices();
+ return se->getMeshElement()->getNumShapeFunctions();
}
virtual int sizeOfC(SElement *se) const
{
- return se->getMeshElement()->getNumVertices();
+ return se->getMeshElement()->getNumShapeFunctions();
}
Dof getLocalDofR(SElement *se, int iRow) const
{
- return Dof(se->getMeshElement()->getVertex(iRow)->getNum(),
+ return Dof(se->getMeshElement()->getShapeFunctionNode(iRow)->getNum(),
Dof::createTypeWithTwoInts(0, _iField));
}
virtual void elementMatrix(SElement *se, fullMatrix<double> &m) const
@@ -41,13 +41,13 @@ class convexCombinationTerm : public femTerm<double> {
MElement *e = se->getMeshElement();
m.setAll(0.);
- const int nbNodes = e->getNumVertices();
+ const int nbSF = e->getNumShapeFunctions();
const double _diff = 1.0;
- for (int j = 0; j < nbNodes; j++){
- for (int k = 0; k < nbNodes; k++){
+ for (int j = 0; j < nbSF; j++){
+ for (int k = 0; k < nbSF; k++){
m(j,k) = -1.*_diff;
}
- m(j,j) = (nbNodes - 1) * _diff;
+ m(j,j) = (nbSF - 1) * _diff;
}
}
@@ -57,11 +57,11 @@ class convexCombinationTerm : public femTerm<double> {
/* MElement *e = se->getMeshElement(); */
/* m.setAll(0.); */
-/* const int nbNodes = e->getNumVertices(); */
-/* for (int j = 0; j < nbNodes; j++){ */
-/* for (int k = 0; k < nbNodes; k++) m(j,k) = 0.0; */
-/* MVertex *v = e->getVertex(j); */
-/* if( v->onWhat()->dim() < 2) m(j,j) = (nbNodes - 1) ; */
+/* const int nbSF = e->getNumShapeFunctions(); */
+/* for (int j = 0; j < nbSF; j++){ */
+/* for (int k = 0; k < nbSF; k++) m(j,k) = 0.0; */
+/* MVertex *v = e->getShapeFunctionNode(j); */
+/* if( v->onWhat()->dim() < 2) m(j,j) = (nbSF - 1) ; */
/* else m(j,j) = 0.0; */
/* } */
/* } */
diff --git a/Solver/crossConfTerm.h b/Solver/crossConfTerm.h
index 011e5306fd..faffd4c086 100644
--- a/Solver/crossConfTerm.h
+++ b/Solver/crossConfTerm.h
@@ -29,26 +29,26 @@ class crossConfTerm : public femTerm<double> {
virtual int sizeOfR(SElement *se) const
{
- return se->getMeshElement()->getNumVertices();
+ return se->getMeshElement()->getNumShapeFunctions();
}
virtual int sizeOfC(SElement *se) const
{
- return se->getMeshElement()->getNumVertices();
+ return se->getMeshElement()->getNumShapeFunctions();
}
Dof getLocalDofR(SElement *se, int iRow) const
{
- return Dof(se->getMeshElement()->getVertex(iRow)->getNum(),
+ return Dof(se->getMeshElement()->getShapeFunctionNode(iRow)->getNum(),
Dof::createTypeWithTwoInts(0, _iFieldR));
}
Dof getLocalDofC(SElement *se, int iRow) const
{
- return Dof(se->getMeshElement()->getVertex(iRow)->getNum(),
+ return Dof(se->getMeshElement()->getShapeFunctionNode(iRow)->getNum(),
Dof::createTypeWithTwoInts(0, _iFieldC));
}
virtual void elementMatrix(SElement *se, fullMatrix<double> &m) const
{
MElement *e = se->getMeshElement();
- int nbNodes = e->getNumVertices();
+ int nbSF = e->getNumShapeFunctions();
int integrationOrder = 2 * (e->getPolynomialOrder() - 1);
int npts;
IntPt *GP;
@@ -70,7 +70,7 @@ class crossConfTerm : public femTerm<double> {
const double _diff = (*_diffusivity)(p.x(), p.y(), p.z());
inv3x3(jac, invjac);
e->getGradShapeFunctions(u, v, w, grads);
- for (int j = 0; j < nbNodes; j++){
+ for (int j = 0; j < nbSF; j++){
Grads[j] = SVector3(invjac[0][0] * grads[j][0] + invjac[0][1] * grads[j][1] +
invjac[0][2] * grads[j][2],
invjac[1][0] * grads[j][0] + invjac[1][1] * grads[j][1] +
@@ -79,42 +79,39 @@ class crossConfTerm : public femTerm<double> {
invjac[2][2] * grads[j][2]);
}
SVector3 N (jac[2][0], jac[2][1], jac[2][2]);
- for (int j = 0; j < nbNodes; j++){
+ for (int j = 0; j < nbSF; j++){
for (int k = 0; k <= j; k++){
m(j, k) += dot(crossprod(Grads[j], Grads[k]), N) * weight * detJ * _diff;
}
}
}
- for (int j = 0; j < nbNodes; j++)
+ for (int j = 0; j < nbSF; j++)
for (int k = 0; k < j; k++)
m(k, j) = -1.* m(j, k);
}
void elementVector(SElement *se, fullVector<double> &m) const
{
-
MElement *e = se->getMeshElement();
- int nbNodes = e->getNumVertices();
+ int nbSF = e->getNumShapeFunctions();
fullMatrix<double> *mat;
- mat = new fullMatrix<double>(nbNodes,nbNodes);
+ mat = new fullMatrix<double>(nbSF,nbSF);
elementMatrix(se, *mat);
- fullVector<double> val(nbNodes);
+ fullVector<double> val(nbSF);
val.scale(0.);
- for (int i = 0; i < nbNodes; i++){
- std::map<MVertex*, SPoint3>::iterator it = _coordView->find(e->getVertex(i));
+ for (int i = 0; i < nbSF; i++){
+ std::map<MVertex*, SPoint3>::iterator it = _coordView->find(e->getShapeFunctionNode(i));
SPoint3 UV = it->second;
if (_iFieldC == 1) val(i) = UV.x();
else if (_iFieldC == 2) val(i) = UV.y();
}
m.scale(0.);
- for (int i = 0; i < nbNodes; i++)
- for (int j = 0; j < nbNodes; j++)
+ for (int i = 0; i < nbSF; i++)
+ for (int j = 0; j < nbSF; j++)
m(i) += -(*mat)(i,j)*val(j);
-
-
}
};
diff --git a/Solver/diagBCTerm.h b/Solver/diagBCTerm.h
index 9b39189829..b7bc918e48 100644
--- a/Solver/diagBCTerm.h
+++ b/Solver/diagBCTerm.h
@@ -24,36 +24,33 @@ class diagBCTerm : public femTerm<double> {
: femTerm<double>(gm), _k(k), _iField(iField) {}
virtual int sizeOfR(SElement *se) const
{
- return se->getMeshElement()->getNumVertices();
+ return se->getMeshElement()->getNumShapeFunctions();
}
virtual int sizeOfC(SElement *se) const
{
- return se->getMeshElement()->getNumVertices();
+ return se->getMeshElement()->getNumShapeFunctions();
}
Dof getLocalDofR(SElement *se, int iRow) const
{
- return Dof(se->getMeshElement()->getVertex(iRow)->getNum(),
+ return Dof(se->getMeshElement()->getShapeFunctionNode(iRow)->getNum(),
Dof::createTypeWithTwoInts(0, _iField));
}
virtual void elementMatrix(SElement *se, fullMatrix<double> &m) const
{
-
MElement *e = se->getMeshElement();
m.setAll(0.);
- const int nbNodes = e->getNumVertices();
- for (int j = 0; j < nbNodes; j++){
- for (int k = 0; k < nbNodes; k++) {
+ const int nbSF = e->getNumShapeFunctions();
+ for (int j = 0; j < nbSF; j++){
+ for (int k = 0; k < nbSF; k++) {
m(j,k) = 0.0;
m(k,j) = 0.0;
}
- MVertex *v = e->getVertex(j);
+ MVertex *v = e->getShapeFunctionNode(j);
if( v->onWhat()->dim() < 2 ) m(j,j) = 1.0;
else m(j,j) = 0.0;
}
- }
-
-
+ }
};
#endif
diff --git a/Solver/distanceTerm.h b/Solver/distanceTerm.h
index dcea0b969b..b8f523e3f5 100644
--- a/Solver/distanceTerm.h
+++ b/Solver/distanceTerm.h
@@ -16,7 +16,6 @@ class distanceTerm : public helmholtzTerm<double> {
void elementVector(SElement *se, fullVector<double> &m) const
{
MElement *e = se->getMeshElement();
- int nbNodes = e->getNumVertices();
int integrationOrder = 2 * e->getPolynomialOrder();
int npts;
IntPt *GP;
@@ -31,7 +30,7 @@ class distanceTerm : public helmholtzTerm<double> {
const double weight = GP[i].weight;
const double detJ = e->getJacobian(u, v, w, jac);
e->getShapeFunctions(u, v, w, ff);
- for (int j = 0; j < nbNodes; j++){
+ for (int j = 0; j < e->getNumShapeFunctions(); j++){
m(j) += ff[j] * weight * detJ;
}
}
diff --git a/Solver/elasticityTerm.cpp b/Solver/elasticityTerm.cpp
index 77a0c4dcca..7c67875577 100644
--- a/Solver/elasticityTerm.cpp
+++ b/Solver/elasticityTerm.cpp
@@ -12,7 +12,7 @@
void elasticityTerm::elementMatrix(SElement *se, fullMatrix<double> &m) const
{
MElement *e = se->getMeshElement();
- int nbNodes = se->getNumNodalShapeFunctions();
+ int nbSF = e->getNumShapeFunctions();
int integrationOrder = 3 * (e->getPolynomialOrder() - 1) ;
int npts;
IntPt *GP;
@@ -43,9 +43,9 @@ void elasticityTerm::elementMatrix(SElement *se, fullMatrix<double> &m) const
{ 0, 0, 0, 0, 0, C44} };
fullMatrix<double> H(6, 6);
- fullMatrix<double> B(6, 3 * nbNodes);
- fullMatrix<double> BTH(3 * nbNodes, 6);
- fullMatrix<double> BT(3 * nbNodes, 6);
+ fullMatrix<double> B(6, 3 * nbSF);
+ fullMatrix<double> BTH(3 * nbSF, 6);
+ fullMatrix<double> BT(3 * nbSF, 6);
for (int i = 0; i < 6; i++)
for (int j = 0; j < 6; j++)
H(i, j) = C[i][j];
@@ -64,35 +64,35 @@ void elasticityTerm::elementMatrix(SElement *se, fullMatrix<double> &m) const
if (se->getShapeEnrichement() == se->getTestEnrichement()){
// printf("coucou\n");
- for (int j = 0; j < nbNodes; j++){
+ for (int j = 0; j < nbSF; j++){
BT(j, 0) = B(0, j) = Grads[j][0];
BT(j, 3) = B(3, j) = Grads[j][1];
BT(j, 5) = B(5, j) = Grads[j][2];
- BT(j + nbNodes, 1) = B(1, j + nbNodes) = Grads[j][1];
- BT(j + nbNodes, 3) = B(3, j + nbNodes) = Grads[j][0];
- BT(j + nbNodes, 4) = B(4, j + nbNodes) = Grads[j][2];
+ BT(j + nbSF, 1) = B(1, j + nbSF) = Grads[j][1];
+ BT(j + nbSF, 3) = B(3, j + nbSF) = Grads[j][0];
+ BT(j + nbSF, 4) = B(4, j + nbSF) = Grads[j][2];
- BT(j + 2 * nbNodes, 2) = B(2, j + 2 * nbNodes) = Grads[j][2];
- BT(j + 2 * nbNodes, 4) = B(4, j + 2 * nbNodes) = Grads[j][1];
- BT(j + 2 * nbNodes, 5) = B(5, j + 2 * nbNodes) = Grads[j][0];
+ BT(j + 2 * nbSF, 2) = B(2, j + 2 * nbSF) = Grads[j][2];
+ BT(j + 2 * nbSF, 4) = B(4, j + 2 * nbSF) = Grads[j][1];
+ BT(j + 2 * nbSF, 5) = B(5, j + 2 * nbSF) = Grads[j][0];
}
}
else{
printf("coucou AAAAAAAAAAAAARGH \n");
se->gradNodalTestFunctions (u, v, w, invjac,GradsT);
- for (int j = 0; j < nbNodes; j++){
+ for (int j = 0; j < nbSF; j++){
BT(j, 0) = Grads[j][0]; B(0, j) = GradsT[j][0];
BT(j, 3) = Grads[j][1]; B(3, j) = GradsT[j][1];
BT(j, 5) = Grads[j][2]; B(5, j) = GradsT[j][2];
- BT(j + nbNodes, 1) = Grads[j][1]; B(1, j + nbNodes) = GradsT[j][1];
- BT(j + nbNodes, 3) = Grads[j][0]; B(3, j + nbNodes) = GradsT[j][0];
- BT(j + nbNodes, 4) = Grads[j][2]; B(4, j + nbNodes) = GradsT[j][2];
+ BT(j + nbSF, 1) = Grads[j][1]; B(1, j + nbSF) = GradsT[j][1];
+ BT(j + nbSF, 3) = Grads[j][0]; B(3, j + nbSF) = GradsT[j][0];
+ BT(j + nbSF, 4) = Grads[j][2]; B(4, j + nbSF) = GradsT[j][2];
- BT(j + 2 * nbNodes, 2) = Grads[j][2]; B(2, j + 2 * nbNodes) = GradsT[j][2];
- BT(j + 2 * nbNodes, 4) = Grads[j][1]; B(4, j + 2 * nbNodes) = GradsT[j][1];
- BT(j + 2 * nbNodes, 5) = Grads[j][0]; B(5, j + 2 * nbNodes) = GradsT[j][0];
+ BT(j + 2 * nbSF, 2) = Grads[j][2]; B(2, j + 2 * nbSF) = GradsT[j][2];
+ BT(j + 2 * nbSF, 4) = Grads[j][1]; B(4, j + 2 * nbSF) = GradsT[j][1];
+ BT(j + 2 * nbSF, 5) = Grads[j][0]; B(5, j + 2 * nbSF) = GradsT[j][0];
}
}
@@ -106,7 +106,7 @@ void elasticityTerm::elementMatrix(SElement *se, fullMatrix<double> &m) const
void elasticityTerm::elementVector(SElement *se, fullVector<double> &m) const
{
MElement *e = se->getMeshElement();
- int nbNodes = e->getNumVertices();
+ int nbSF = e->getNumShapeFunctions();
int integrationOrder = 2 * e->getPolynomialOrder();
int npts;
IntPt *GP;
@@ -123,10 +123,10 @@ void elasticityTerm::elementVector(SElement *se, fullVector<double> &m) const
const double weight = GP[i].weight;
const double detJ = e->getJacobian(u, v, w, jac);
se->nodalTestFunctions(u, v, w, ff);
- for (int j = 0; j < nbNodes; j++){
+ for (int j = 0; j < nbSF; j++){
m(j) += _volumeForce.x() * ff[j] * weight * detJ * .5;
- m(j + nbNodes) += _volumeForce.y() * ff[j] * weight * detJ * .5;
- m(j + 2 * nbNodes) += _volumeForce.z() * ff[j] * weight * detJ * .5;
+ m(j + nbSF) += _volumeForce.y() * ff[j] * weight * detJ * .5;
+ m(j + 2 * nbSF) += _volumeForce.z() * ff[j] * weight * detJ * .5;
}
}
}
diff --git a/Solver/elasticityTerm.h b/Solver/elasticityTerm.h
index 9d131a01b0..c89b1593d7 100644
--- a/Solver/elasticityTerm.h
+++ b/Solver/elasticityTerm.h
@@ -23,28 +23,28 @@ class elasticityTerm : public femTerm<double> {
// element matrix size : 3 dofs per vertex
virtual int sizeOfR(SElement *se) const
{
- return 3 * se->getMeshElement()->getNumVertices();
+ return 3 * se->getMeshElement()->getNumShapeFunctions();
}
virtual int sizeOfC(SElement *se) const
{
- return 3 * se->getMeshElement()->getNumVertices();
+ return 3 * se->getMeshElement()->getNumShapeFunctions();
}
// order dofs in the local matrix :
// dx1, dx2 ... dxn, dy1, ..., dyn, dz1, ... dzn
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(),
+ int iCompR = iRow / e->getNumShapeFunctions();
+ int ithLocalVertex = iRow % e->getNumShapeFunctions();
+ return Dof(e->getShapeFunctionNode(ithLocalVertex)->getNum(),
Dof::createTypeWithTwoInts(iCompR, _iFieldR));
}
Dof getLocalDofC(SElement *se, int iCol) const
{
MElement *e = se->getMeshElement();
- int iCompC = iCol / e->getNumVertices();
- int ithLocalVertex = iCol % e->getNumVertices();
- return Dof(e->getVertex(ithLocalVertex)->getNum(),
+ int iCompC = iCol / e->getNumShapeFunctions();
+ int ithLocalVertex = iCol % e->getNumShapeFunctions();
+ return Dof(e->getShapeFunctionNode(ithLocalVertex)->getNum(),
Dof::createTypeWithTwoInts(iCompC, _iFieldC));
}
public:
diff --git a/Solver/functionSpace.cpp b/Solver/functionSpace.cpp
index 6cb8e84b28..09fc017e23 100644
--- a/Solver/functionSpace.cpp
+++ b/Solver/functionSpace.cpp
@@ -11,5 +11,7 @@
//
#include "functionSpace.h"
-const SVector3 VectorLagrangeFunctionSpaceOfElement::BasisVectors[3]={SVector3(1,0,0),SVector3(0,1,0),SVector3(0,0,1)};
-const SVector3 VectorLagrangeFunctionSpace::BasisVectors[3]={SVector3(1,0,0),SVector3(0,1,0),SVector3(0,0,1)};
+const SVector3 VectorLagrangeFunctionSpaceOfElement::BasisVectors[3] =
+ {SVector3(1, 0, 0), SVector3(0, 1, 0), SVector3(0, 0, 1)};
+const SVector3 VectorLagrangeFunctionSpace::BasisVectors[3] =
+ {SVector3(1, 0, 0), SVector3(0, 1, 0), SVector3(0, 0, 1)};
diff --git a/Solver/functionSpace.h b/Solver/functionSpace.h
index 1568c82095..644378f027 100644
--- a/Solver/functionSpace.h
+++ b/Solver/functionSpace.h
@@ -54,7 +54,7 @@ class FunctionSpace : public FunctionSpaceBase
typedef typename TensorialTraits<T>::GradType GradType;
typedef typename TensorialTraits<T>::HessType HessType;
virtual void f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals) = 0;
- virtual void fuvw(MElement *ele, double u, double v, double w, std::vector<ValType> &vals) {}; // should return to pure virtual once all is done.
+ virtual void fuvw(MElement *ele, double u, double v, double w, std::vector<ValType> &vals) {} // should return to pure virtual once all is done.
virtual void gradf(MElement *ele, double u, double v, double w, std::vector<GradType> &grads) = 0;
virtual void gradfuvw(MElement *ele, double u, double v, double w, std::vector<GradType> &grads) {} // should return to pure virtual once all is done.
virtual void hessfuvw(MElement *ele, double u, double v, double w, std::vector<HessType> &hess) = 0;
@@ -88,7 +88,7 @@ class ScalarLagrangeFunctionSpaceOfElement : public FunctionSpace<double>
ele->movePointFromParentSpaceToElementSpace(u, v, w);
}
}
- int ndofs = ele->getNumVertices();
+ int ndofs = ele->getNumShapeFunctions();
int curpos = vals.size();
vals.resize(curpos + ndofs);
ele->getShapeFunctions(u, v, w, &(vals[curpos]));
@@ -101,7 +101,7 @@ class ScalarLagrangeFunctionSpaceOfElement : public FunctionSpace<double>
ele->movePointFromParentSpaceToElementSpace(u, v, w);
}
}
- int ndofs = ele->getNumVertices();
+ int ndofs = ele->getNumShapeFunctions();
grads.reserve(grads.size() + ndofs);
double gradsuvw[256][3];
ele->getGradShapeFunctions(u, v, w, gradsuvw);
@@ -109,7 +109,7 @@ class ScalarLagrangeFunctionSpaceOfElement : public FunctionSpace<double>
double invjac[3][3];
const double detJ = ele->getJacobian(u, v, w, jac); // redondant : on fait cet appel a l'exterieur
inv3x3(jac, invjac);
- for (int i = 0; i < ndofs; ++i)
+ for(int i = 0; i < ndofs; ++i)
grads.push_back(GradType(
invjac[0][0] * gradsuvw[i][0] + invjac[0][1] * gradsuvw[i][1] + invjac[0][2] * gradsuvw[i][2],
invjac[1][0] * gradsuvw[i][0] + invjac[1][1] * gradsuvw[i][1] + invjac[1][2] * gradsuvw[i][2],
@@ -123,19 +123,18 @@ class ScalarLagrangeFunctionSpaceOfElement : public FunctionSpace<double>
ele->movePointFromParentSpaceToElementSpace(u, v, w);
}
}
- int ndofs = ele->getNumVertices(); // ATTENTION RETOURNE LE NBBRE DE NOEUDS ET PAS LE NBRE DE DDL
+ int ndofs = ele->getNumShapeFunctions();
hess.reserve(hess.size() + ndofs); // permet de mettre les composantes suivantes à la suite du vecteur
double hessuvw[256][3][3];
ele->getHessShapeFunctions(u, v, w, hessuvw);
HessType hesst;
- for (int i = 0; i < ndofs; ++i){
+ for(int i = 0; i < ndofs; ++i){
hesst(0,0) = hessuvw[i][0][0]; hesst(0,1) = hessuvw[i][0][1]; hesst(0,2) = hessuvw[i][0][2];
hesst(1,0) = hessuvw[i][1][0]; hesst(1,1) = hessuvw[i][1][1]; hesst(1,2) = hessuvw[i][1][2];
hesst(2,0) = hessuvw[i][2][0]; hesst(2,1) = hessuvw[i][2][1]; hesst(2,2) = hessuvw[i][2][2];
hess.push_back(hesst);
}
}
-
virtual void gradfuvw(MElement *ele, double u, double v, double w, std::vector<GradType> &grads)
{
if(ele->getParent()) {
@@ -143,25 +142,23 @@ class ScalarLagrangeFunctionSpaceOfElement : public FunctionSpace<double>
ele->movePointFromParentSpaceToElementSpace(u, v, w);
}
}
- int ndofs = ele->getNumVertices();
+ int ndofs = ele->getNumShapeFunctions();
grads.reserve(grads.size() + ndofs);
double gradsuvw[256][3];
ele->getGradShapeFunctions(u, v, w, gradsuvw);
- for (int i = 0; i < ndofs; ++i)
+ for(int i = 0; i < ndofs; ++i)
grads.push_back(GradType(gradsuvw[i][0], gradsuvw[i][1], gradsuvw[i][2]));
}
-
virtual int getNumKeys(MElement *ele)
{
- return ele->getNumVertices();
+ return ele->getNumShapeFunctions();
}
-
virtual void getKeys(MElement *ele, std::vector<Dof> &keys) // appends ...
{
- int ndofs = ele->getNumVertices();
+ int ndofs = ele->getNumShapeFunctions();
keys.reserve(keys.size() + ndofs);
- for (int i = 0; i < ndofs; ++i)
- getKeys(ele->getVertex(i), keys);
+ for(int i = 0; i < ndofs; ++i)
+ getKeys(ele->getShapeFunctionNode(i), keys);
}
};
@@ -186,17 +183,16 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
virtual void f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)
{
if(ele->getParent()) ele = ele->getParent();
- int ndofs = ele->getNumVertices();
+ int ndofs = ele->getNumShapeFunctions();
int curpos = vals.size();
vals.resize(curpos + ndofs);
ele->getShapeFunctions(u, v, w, &(vals[curpos]));
}
-
// Fonction renvoyant un vecteur contenant le grandient de chaque FF
virtual void gradf(MElement *ele, double u, double v, double w, std::vector<GradType> &grads)
{
if(ele->getParent()) ele = ele->getParent();
- int ndofs = ele->getNumVertices();
+ int ndofs = ele->getNumShapeFunctions();
grads.reserve(grads.size() + ndofs);
double gradsuvw[256][3];
ele->getGradShapeFunctions(u, v, w, gradsuvw);
@@ -204,22 +200,22 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
double invjac[3][3];
const double detJ = ele->getJacobian(u, v, w, jac); // redondant : on fait cet appel a l'exterieur
inv3x3(jac, invjac);
- for (int i = 0; i < ndofs; ++i)
+ for(int i = 0; i < ndofs; ++i)
grads.push_back(GradType(
invjac[0][0] * gradsuvw[i][0] + invjac[0][1] * gradsuvw[i][1] + invjac[0][2] * gradsuvw[i][2],
invjac[1][0] * gradsuvw[i][0] + invjac[1][1] * gradsuvw[i][1] + invjac[1][2] * gradsuvw[i][2],
invjac[2][0] * gradsuvw[i][0] + invjac[2][1] * gradsuvw[i][1] + invjac[2][2] * gradsuvw[i][2]));
}
- // Fonction renvoyant un vecteur contenant le hessien [][] de chaque FF dans l'espace ISOPARAMETRIQUE
+ // Fonction renvoyant un vecteur contenant le hessien [][] de chaque FF dans l'espace ISOPARAMETRIQUE
virtual void hessfuvw(MElement *ele, double u, double v, double w, std::vector<HessType> &hess)
{
if(ele->getParent()) ele = ele->getParent();
- int ndofs = ele->getNumVertices(); // ATTENTION RETOURNE LE NBBRE DE NOEUDS ET PAS LE NBRE DE DDL
+ int ndofs = ele->getNumShapeFunctions();
hess.reserve(hess.size() + ndofs); // permet de mettre les composantes suivantes à la suite du vecteur
double hessuvw[256][3][3];
ele->getHessShapeFunctions(u, v, w, hessuvw);
HessType hesst;
- for (int i = 0; i < ndofs; ++i){
+ for(int i = 0; i < ndofs; ++i){
hesst(0,0) = hessuvw[i][0][0]; hesst(0,1) = hessuvw[i][0][1]; hesst(0,2) = hessuvw[i][0][2];
hesst(1,0) = hessuvw[i][1][0]; hesst(1,1) = hessuvw[i][1][1]; hesst(1,2) = hessuvw[i][1][2];
hesst(2,0) = hessuvw[i][2][0]; hesst(2,1) = hessuvw[i][2][1]; hesst(2,2) = hessuvw[i][2][2];
@@ -229,37 +225,35 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
virtual void gradfuvw(MElement *ele, double u, double v, double w, std::vector<GradType> &grads)
{
if(ele->getParent()) ele = ele->getParent();
- int ndofs = ele->getNumVertices();
+ int ndofs = ele->getNumShapeFunctions();
grads.reserve(grads.size() + ndofs);
double gradsuvw[256][3];
ele->getGradShapeFunctions(u, v, w, gradsuvw);
- for (int i = 0; i < ndofs; ++i)
+ for(int i = 0; i < ndofs; ++i)
grads.push_back(GradType(gradsuvw[i][0], gradsuvw[i][1], gradsuvw[i][2]));
}
-
virtual void fuvw(MElement *ele, double u, double v, double w,std::vector<ValType> &vals)
{
- if (ele->getParent()) ele = ele->getParent();
- int ndofs= ele->getNumVertices();
+ if(ele->getParent()) ele = ele->getParent();
+ int ndofs= ele->getNumShapeFunctions();
vals.reserve(vals.size()+ndofs);
double valsuvw[256];
ele->getShapeFunctions(u, v, w, valsuvw);
- for (int i=0;i<ndofs;++i) {vals.push_back(valsuvw[i]);}
+ for(int i = 0; i < ndofs; ++i)
+ vals.push_back(valsuvw[i]);
}
-
virtual int getNumKeys(MElement *ele)
{
if(ele->getParent()) ele = ele->getParent();
- return ele->getNumVertices();
+ return ele->getNumShapeFunctions();
}
-
virtual void getKeys(MElement *ele, std::vector<Dof> &keys) // appends ...
{
if(ele->getParent()) ele = ele->getParent();
- int ndofs = ele->getNumVertices();
+ int ndofs = ele->getNumShapeFunctions();
keys.reserve(keys.size() + ndofs);
- for (int i = 0; i < ndofs; ++i)
- getKeys(ele->getVertex(i), keys);
+ for(int i = 0; i < ndofs; ++i)
+ getKeys(ele->getShapeFunctionNode(i), keys);
}
};
@@ -283,7 +277,7 @@ public :
const T& mult2, int comp2_): ScalarFS(new T2(SFS))
{
multipliers.push_back(mult1); multipliers.push_back(mult2);
- comp.push_back(comp1_); comp.push_back(comp2_);
+ comp.push_back(comp1_); comp.push_back(comp2_);
}
template <class T2> ScalarToAnyFunctionSpace(const T2 &SFS, const T& mult1, int comp1_,
@@ -303,9 +297,8 @@ public :
int nbcomp = comp.size();
int curpos = vals.size();
vals.reserve(curpos + nbcomp * nbdofs);
- for (int j = 0; j < nbcomp; ++j)
- {
- for (int i = 0; i < nbdofs; ++i)
+ for(int j = 0; j < nbcomp; ++j){
+ for(int i = 0; i < nbdofs; ++i)
vals.push_back(multipliers[j] * valsd[i]);
}
}
@@ -319,10 +312,8 @@ public :
int curpos = grads.size();
grads.reserve(curpos + nbcomp * nbdofs);
GradType val;
- for (int j = 0; j < nbcomp; ++j)
- {
- for (int i = 0; i < nbdofs; ++i)
- {
+ for(int j = 0; j < nbcomp; ++j){
+ for(int i = 0; i < nbdofs; ++i){
tensprod(multipliers[j], gradsd[i], val);
grads.push_back(val);
}
@@ -341,10 +332,8 @@ public :
int curpos = grads.size();
grads.reserve(curpos + nbcomp * nbdofs);
GradType val;
- for (int j = 0; j < nbcomp; ++j)
- {
- for (int i = 0; i < nbdofs; ++i)
- {
+ for(int j = 0; j < nbcomp; ++j){
+ for(int i = 0; i < nbdofs; ++i){
tensprod(multipliers[j], gradsd[i], val);
grads.push_back(val);
}
@@ -363,10 +352,8 @@ public :
int nbcomp = comp.size();
int curpos = keys.size();
keys.reserve(curpos + nbcomp * nbdofs);
- for (int j = 0; j < nbcomp; ++j)
- {
- for (int i = 0; i < nk; ++i)
- {
+ for(int j = 0; j < nbcomp; ++j){
+ for(int i = 0; i < nk; ++i){
int i1, i2;
Dof::getTwoIntsFromType(bufk[i].getType(), i1, i2);
keys.push_back(Dof(bufk[i].getEntity(), Dof::createTypeWithTwoInts(comp[j], i1)));
@@ -387,15 +374,15 @@ class VectorLagrangeFunctionSpaceOfElement : public ScalarToAnyFunctionSpace<SVe
ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeFunctionSpaceOfElement(id),
SVector3(1.,0.,0.), VECTOR_X, SVector3(0.,1.,0.), VECTOR_Y, SVector3(0.,0.,1.), VECTOR_Z)
{}
- VectorLagrangeFunctionSpaceOfElement(int id,Along comp1) :
+ VectorLagrangeFunctionSpaceOfElement(int id, Along comp1) :
ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeFunctionSpaceOfElement(id),
BasisVectors[comp1], comp1)
{}
- VectorLagrangeFunctionSpaceOfElement(int id,Along comp1,Along comp2) :
+ VectorLagrangeFunctionSpaceOfElement(int id, Along comp1, Along comp2) :
ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeFunctionSpaceOfElement(id),
BasisVectors[comp1], comp1, BasisVectors[comp2], comp2)
{}
- VectorLagrangeFunctionSpaceOfElement(int id,Along comp1,Along comp2, Along comp3) :
+ VectorLagrangeFunctionSpaceOfElement(int id, Along comp1, Along comp2, Along comp3) :
ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeFunctionSpaceOfElement(id),
BasisVectors[comp1], comp1, BasisVectors[comp2], comp2, BasisVectors[comp3], comp3)
{}
@@ -413,15 +400,15 @@ class VectorLagrangeFunctionSpace : public ScalarToAnyFunctionSpace<SVector3>
ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeFunctionSpace(id),
SVector3(1.,0.,0.), VECTOR_X, SVector3(0.,1.,0.), VECTOR_Y, SVector3(0.,0.,1.), VECTOR_Z)
{}
- VectorLagrangeFunctionSpace(int id,Along comp1) :
+ VectorLagrangeFunctionSpace(int id, Along comp1) :
ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeFunctionSpace(id),
BasisVectors[comp1], comp1)
{}
- VectorLagrangeFunctionSpace(int id,Along comp1,Along comp2) :
+ VectorLagrangeFunctionSpace(int id, Along comp1, Along comp2) :
ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeFunctionSpace(id),
BasisVectors[comp1], comp1, BasisVectors[comp2], comp2)
{}
- VectorLagrangeFunctionSpace(int id,Along comp1,Along comp2, Along comp3) :
+ VectorLagrangeFunctionSpace(int id, Along comp1, Along comp2, Along comp3) :
ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeFunctionSpace(id),
BasisVectors[comp1], comp1, BasisVectors[comp2], comp2, BasisVectors[comp3], comp3)
{}
@@ -440,7 +427,7 @@ class CompositeFunctionSpace : public FunctionSpace<T>
std::vector<FunctionSpace<T>* > _spaces;
public:
template <class T1> CompositeFunctionSpace(const T1& t) { _spaces.push_back(new T1(t));}
- template <class T1, class T2> CompositeFunctionSpace(const T1& t1,const T2& t2)
+ template <class T1, class T2> CompositeFunctionSpace(const T1& t1, const T2& t2)
{ _spaces.push_back(new T1(t1));
_spaces.push_back(new T2(t2)); }
template <class T1, class T2, class T3> CompositeFunctionSpace(const T1& t1, const T2& t2, const T3& t3)
@@ -527,17 +514,16 @@ class FilteredFunctionSpace : public FunctionSpace<T>
FunctionSpace<T>* _spacebase;
F *_filter;
public:
- virtual void hessfuvw(MElement *ele, double u, double v, double w,std::vector<HessType> &hess){}
- FilteredFunctionSpace<T,F>(FunctionSpace<T>* spacebase,F * filter) : _spacebase(spacebase),_filter(filter){}
- virtual void f(MElement *ele, double u, double v, double w,std::vector<ValType> &vals);
- virtual void gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads);
+ virtual void hessfuvw(MElement *ele, double u, double v, double w, std::vector<HessType> &hess){}
+ FilteredFunctionSpace<T,F>(FunctionSpace<T>* spacebase,F * filter) : _spacebase(spacebase), _filter(filter){}
+ virtual void f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals);
+ virtual void gradf(MElement *ele, double u, double v, double w, std::vector<GradType> &grads);
virtual int getNumKeys(MElement *ele);
virtual void getKeys(MElement *ele, std::vector<Dof> &keys);
};
-
-template <class T> void xFemFunctionSpace<T>::f(MElement *ele, double u, double v, double w,std::vector<ValType> &vals)
+template <class T> void xFemFunctionSpace<T>::f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)
{
// We need parent parameters
MElement * elep;
@@ -546,78 +532,72 @@ template <class T> void xFemFunctionSpace<T>::f(MElement *ele, double u, double
// Get the spacebase valsd
std::vector<ValType> valsd;
- xFemFunctionSpace<T>::_spacebase->f(elep,u,v,w,valsd);
+ xFemFunctionSpace<T>::_spacebase->f(elep, u, v, w, valsd);
- int nbdofs=valsd.size();
- int curpos=vals.size(); // if in composite function space
- vals.reserve(curpos+nbdofs);
+ int nbdofs = valsd.size();
+ int curpos = vals.size(); // if in composite function space
+ vals.reserve(curpos + nbdofs);
- // then enriched dofs so the order is ex:(a2x,a2y,a3x,a3y)
- if (nbdofs>0) // if enriched
- {
- // Enrichment function calcul
+ // then enriched dofs so the order is ex:(a2x,a2y,a3x,a3y)
+ if (nbdofs > 0){ // if enriched
+ // Enrichment function calcul
SPoint3 p;
elep->pnt(u, v, w, p); // parametric to cartesian coordinates
double func;
_funcEnrichment->setElement(elep);
func = (*_funcEnrichment)(p.x(), p.y(), p.z());
- for (int i=0 ;i<nbdofs;i++)
- {
- vals.push_back(valsd[i]*func);
+ for (int i = 0 ; i < nbdofs; i++){
+ vals.push_back(valsd[i] * func);
}
}
}
-template <class T> void xFemFunctionSpace<T>::gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
+template <class T> void xFemFunctionSpace<T>::gradf(MElement *ele, double u, double v, double w, std::vector<GradType> &grads)
{
// We need parent parameters
- MElement * elep;
+ MElement *elep;
if (ele->getParent()) elep = ele->getParent();
else elep = ele;
// Get the spacebase gradsd
std::vector<GradType> gradsd;
- xFemFunctionSpace<T>::_spacebase->gradf(elep,u,v,w,gradsd);
-
+ xFemFunctionSpace<T>::_spacebase->gradf(elep, u, v, w, gradsd);
int nbdofs=gradsd.size();
// We need spacebase valsd to compute total gradient
std::vector<ValType> valsd;
- xFemFunctionSpace<T>::_spacebase->f(elep,u,v,w,valsd);
+ xFemFunctionSpace<T>::_spacebase->f(elep, u, v, w, valsd);
- int curpos=grads.size(); // if in composite function space
- grads.reserve(curpos+nbdofs);
+ int curpos = grads.size(); // if in composite function space
+ grads.reserve(curpos + nbdofs);
// then enriched dofs so the order is ex:(a2x,a2y,a3x,a3y)
- if (nbdofs>0) // if enriched
- {
+ if (nbdofs > 0){ // if enriched
double df[3];
SPoint3 p;
elep->pnt(u, v, w, p);
_funcEnrichment->setElement(elep);
- _funcEnrichment->gradient (p.x(), p.y(),p.z(),df[0],df[1],df[2]);
- ValType gradfuncenrich(df[0],df[1],df[2]);
+ _funcEnrichment->gradient (p.x(), p.y(), p.z(), df[0], df[1], df[2]);
+ ValType gradfuncenrich(df[0], df[1], df[2]);
// Enrichment function calcul
-
double func;
_funcEnrichment->setElement(elep);
func = (*_funcEnrichment)(p.x(), p.y(), p.z());
- for (int i=0 ;i<nbdofs;i++)
- {
+ for (int i = 0 ; i < nbdofs; i++){
GradType GradFunc;
- tensprod(valsd[i],gradfuncenrich,GradFunc);
- grads.push_back(gradsd[i]*func + GradFunc);
+ tensprod(valsd[i], gradfuncenrich, GradFunc);
+ grads.push_back(gradsd[i] * func + GradFunc);
}
}
}
template <class T> int xFemFunctionSpace<T>::getNumKeys(MElement *ele)
{
- MElement * elep;
+ MElement *elep;
if (ele->getParent()) elep = ele->getParent();
else elep = ele;
int nbdofs = xFemFunctionSpace<T>::_spacebase->getNumKeys(elep);
@@ -626,28 +606,27 @@ template <class T> int xFemFunctionSpace<T>::getNumKeys(MElement *ele)
template <class T> void xFemFunctionSpace<T>::getKeys(MElement *ele, std::vector<Dof> &keys)
{
- MElement * elep;
+ MElement *elep;
if (ele->getParent()) elep = ele->getParent();
else elep = ele;
- int normalk=xFemFunctionSpace<T>::_spacebase->getNumKeys(elep);
+ int normalk = xFemFunctionSpace<T>::_spacebase->getNumKeys(elep);
std::vector<Dof> bufk;
bufk.reserve(normalk);
- xFemFunctionSpace<T>::_spacebase->getKeys(elep,bufk);
- int normaldofs=bufk.size();
+ xFemFunctionSpace<T>::_spacebase->getKeys(elep, bufk);
+ int normaldofs = bufk.size();
int nbdofs = normaldofs;
- int curpos=keys.size();
- keys.reserve(curpos+nbdofs);
+ int curpos = keys.size();
+ keys.reserve(curpos + nbdofs);
// get keys so the order is ex:(a2x,a2y,a3x,a3y)
// enriched dof tagged with ( i2 -> i2 + 1 )
- for (int i=0 ;i<nbdofs;i++)
- {
- int i1,i2;
- Dof::getTwoIntsFromType(bufk[i].getType(), i1,i2);
- keys.push_back(Dof(bufk[i].getEntity(),Dof::createTypeWithTwoInts(i1,i2+1)));
+ for (int i = 0 ;i < nbdofs; i++){
+ int i1, i2;
+ Dof::getTwoIntsFromType(bufk[i].getType(), i1, i2);
+ keys.push_back(Dof(bufk[i].getEntity(), Dof::createTypeWithTwoInts(i1, i2 + 1)));
}
}
@@ -655,54 +634,50 @@ template <class T> void xFemFunctionSpace<T>::getKeys(MElement *ele, std::vector
// Filtered function space
//
-template <class T,class F> void FilteredFunctionSpace<T,F>::f(MElement *ele, double u, double v, double w,std::vector<ValType> &vals)
+template <class T,class F> void FilteredFunctionSpace<T,F>::f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)
{
// We need parent parameters
- MElement * elep;
+ MElement *elep;
if (ele->getParent()) elep = ele->getParent();
else elep = ele;
std::vector<ValType> valsd;
- _spacebase->f(elep,u,v,w,valsd);
+ _spacebase->f(elep, u, v, w, valsd);
- int normalk=_spacebase->getNumKeys(elep);
+ int normalk = _spacebase->getNumKeys(elep);
std::vector<Dof> bufk;
bufk.reserve(normalk);
- _spacebase->getKeys(elep,bufk);
+ _spacebase->getKeys(elep, bufk);
- for (int i=0;i<bufk.size();i++)
- {
+ for (int i = 0; i < bufk.size(); i++){
if ((*_filter)(bufk[i]))
vals.push_back(valsd[i]);
}
-
}
-template <class T,class F> void FilteredFunctionSpace<T,F>::gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
+template <class T,class F> void FilteredFunctionSpace<T,F>::gradf(MElement *ele, double u, double v, double w, std::vector<GradType> &grads)
{
// We need parent parameters
- MElement * elep;
+ MElement *elep;
if (ele->getParent()) elep = ele->getParent();
else elep = ele;
// Get space base gradsd
std::vector<GradType> gradsd;
- _spacebase->gradf(elep,u,v,w,gradsd);
+ _spacebase->gradf(elep, u, v, w, gradsd);
// Get numkeys
- int normalk=_spacebase->getNumKeys(elep);
+ int normalk = _spacebase->getNumKeys(elep);
std::vector<Dof> bufk;
bufk.reserve(normalk);
- _spacebase->getKeys(elep,bufk);
+ _spacebase->getKeys(elep, bufk);
- for (int i=0;i<bufk.size();i++)
- {
- if ( (*_filter)(bufk[i]))
+ for (int i = 0; i < bufk.size(); i++){
+ if ((*_filter)(bufk[i]))
grads.push_back(gradsd[i]);
}
-
}
template <class T,class F> int FilteredFunctionSpace<T,F>::getNumKeys(MElement *ele)
@@ -713,36 +688,34 @@ template <class T,class F> int FilteredFunctionSpace<T,F>::getNumKeys(MElement *
int nbdofs = 0;
- int normalk=_spacebase->getNumKeys(elep);
+ int normalk = _spacebase->getNumKeys(elep);
std::vector<Dof> bufk;
bufk.reserve(normalk);
- _spacebase->getKeys(elep,bufk);
+ _spacebase->getKeys(elep, bufk);
- for (int i=0;i<bufk.size();i++)
- {
- if ( (*_filter)(bufk[i]))
+ for (int i = 0; i < bufk.size(); i++){
+ if ((*_filter)(bufk[i]))
nbdofs = nbdofs + 1;
}
return nbdofs;
}
-template <class T,class F> void FilteredFunctionSpace<T,F>::getKeys(MElement *ele, std::vector<Dof> &keys)
+template <class T, class F> void FilteredFunctionSpace<T, F>::getKeys(MElement *ele, std::vector<Dof> &keys)
{
- MElement * elep;
+ MElement *elep;
if (ele->getParent()) elep = ele->getParent();
else elep = ele;
- int normalk=_spacebase->getNumKeys(elep);
+ int normalk = _spacebase->getNumKeys(elep);
std::vector<Dof> bufk;
bufk.reserve(normalk);
- _spacebase->getKeys(elep,bufk);
- int normaldofs=bufk.size();
+ _spacebase->getKeys(elep, bufk);
+ int normaldofs = bufk.size();
int nbdofs = 0;
- for (int i=0;i<bufk.size();i++)
- {
- if ( (*_filter)(bufk[i]))
+ for (int i = 0; i < bufk.size(); i++){
+ if ((*_filter)(bufk[i]))
keys.push_back(bufk[i]);
}
}
diff --git a/Solver/helmholtzTerm.h b/Solver/helmholtzTerm.h
index 5946d62e68..da781d22f5 100644
--- a/Solver/helmholtzTerm.h
+++ b/Solver/helmholtzTerm.h
@@ -29,20 +29,20 @@ class helmholtzTerm : public femTerm<scalar> {
// one dof per vertex (nodal fem)
virtual int sizeOfR(SElement *se) const
{
- return se->getMeshElement()->getNumVertices();
+ return se->getMeshElement()->getNumShapeFunctions();
}
virtual int sizeOfC(SElement *se) const
{
- return se->getMeshElement()->getNumVertices();
+ return se->getMeshElement()->getNumShapeFunctions();
}
Dof getLocalDofR(SElement *se, int iRow) const
{
- return Dof(se->getMeshElement()->getVertex(iRow)->getNum(),
+ return Dof(se->getMeshElement()->getShapeFunctionNode(iRow)->getNum(),
Dof::createTypeWithTwoInts(0, _iFieldR));
}
Dof getLocalDofC(SElement *se, int iRow) const
{
- return Dof(se->getMeshElement()->getVertex(iRow)->getNum(),
+ return Dof(se->getMeshElement()->getShapeFunctionNode(iRow)->getNum(),
Dof::createTypeWithTwoInts(0, _iFieldC));
}
virtual void elementMatrix(SElement *se, fullMatrix<scalar> &m) const
@@ -57,9 +57,9 @@ class helmholtzTerm : public femTerm<scalar> {
int npts; IntPt *GP;
e->getIntegrationPoints(integrationOrder, &npts, &GP);
// get the number of nodes
- const int nbNodes = e->getNumVertices();
+ const int nbSF = e->getNumShapeFunctions();
// assume a maximum of 100 nodes
- assert(nbNodes < 100);
+ assert(nbSF < 100);
double jac[3][3];
double invjac[3][3];
double Grads[100][3], grads[100][3];
@@ -79,7 +79,7 @@ class helmholtzTerm : public femTerm<scalar> {
inv3x3(jac, invjac) ;
e->getGradShapeFunctions(u, v, w, grads);
if (_a) e->getShapeFunctions(u, v, w, sf);
- for (int j = 0; j < nbNodes; j++){
+ for (int j = 0; j < nbSF; j++){
Grads[j][0] = invjac[0][0] * grads[j][0] + invjac[0][1] * grads[j][1] +
invjac[0][2] * grads[j][2];
Grads[j][1] = invjac[1][0] * grads[j][0] + invjac[1][1] * grads[j][1] +
@@ -88,7 +88,7 @@ class helmholtzTerm : public femTerm<scalar> {
invjac[2][2] * grads[j][2];
if (!_a) sf[j] = 0;
}
- for (int j = 0; j < nbNodes; j++){
+ for (int j = 0; j < nbSF; j++){
for (int k = 0; k <= j; k++){
m(j, k) += (K * (Grads[j][0] * Grads[k][0] +
Grads[j][1] * Grads[k][1] +
@@ -96,7 +96,7 @@ class helmholtzTerm : public femTerm<scalar> {
}
}
}
- for (int j = 0; j < nbNodes; j++)
+ for (int j = 0; j < nbSF; j++)
for (int k = 0; k < j; k++)
m(k, j) = m(j, k);
diff --git a/Solver/laplaceTerm.h b/Solver/laplaceTerm.h
index 7c8f29430b..94f44cf59b 100644
--- a/Solver/laplaceTerm.h
+++ b/Solver/laplaceTerm.h
@@ -19,28 +19,26 @@ class laplaceTerm : public helmholtzTerm<double> {
: helmholtzTerm<double>(gm, iField, iField, k, 0), _iField(iField), _coordView(coord) {}
void elementVector(SElement *se, fullVector<double> &m) const
{
-
MElement *e = se->getMeshElement();
- int nbNodes = e->getNumVertices();
+ int nbSF = e->getNumShapeFunctions();
fullMatrix<double> *mat;
- mat = new fullMatrix<double>(nbNodes,nbNodes);
+ mat = new fullMatrix<double>(nbSF,nbSF);
elementMatrix(se, *mat);
- fullVector<double> val(nbNodes);
+ fullVector<double> val(nbSF);
val.scale(0.);
- for (int i = 0; i < nbNodes; i++){
- std::map<MVertex*, SPoint3>::iterator it = _coordView->find(e->getVertex(i));
+ for (int i = 0; i < nbSF; i++){
+ std::map<MVertex*, SPoint3>::iterator it = _coordView->find(e->getShapeFunctionNode(i));
SPoint3 UV = it->second;
if (_iField == 1) val(i) = UV.x();
else if (_iField == 2) val(i) = UV.y();
}
m.scale(0.);
- for (int i = 0; i < nbNodes; i++)
- for (int j = 0; j < nbNodes; j++)
- m(i) += -(*mat)(i,j)*val(j);
-
+ for (int i = 0; i < nbSF; i++)
+ for (int j = 0; j < nbSF; j++)
+ m(i) += -(*mat)(i,j) * val(j);
}
};
diff --git a/Solver/orthogonalTerm.h b/Solver/orthogonalTerm.h
index 6a5bd80ac9..908f6c624d 100644
--- a/Solver/orthogonalTerm.h
+++ b/Solver/orthogonalTerm.h
@@ -24,12 +24,11 @@ class orthogonalTerm : public helmholtzTerm<double> {
: helmholtzTerm<double>(gm, iField, iField, k, 0), _dofView(dofView), withDof(true) {}
void elementVector(SElement *se, fullVector<double> &m) const
{
-
MElement *e = se->getMeshElement();
- int nbNodes = e->getNumVertices();
+ int nbSF = e->getNumShapeFunctions();
//fill elementary matrix mat(i,j)
- int integrationOrder = 2* (e->getPolynomialOrder() - 1);
+ int integrationOrder = 2 * (e->getPolynomialOrder() - 1);
int npts;
IntPt *GP;
double jac[3][3];
@@ -37,10 +36,10 @@ class orthogonalTerm : public helmholtzTerm<double> {
SVector3 Grads [256];
double grads[256][3];
e->getIntegrationPoints(integrationOrder, &npts, &GP);
- fullMatrix<double> mat(nbNodes,nbNodes);
+ fullMatrix<double> mat(nbSF, nbSF);
mat.setAll(0.);
- for (int i = 0; i < npts; i++){
+ for(int i = 0; i < npts; i++){
const double u = GP[i].pt[0];
const double v = GP[i].pt[1];
const double w = GP[i].pt[2];
@@ -49,7 +48,7 @@ class orthogonalTerm : public helmholtzTerm<double> {
SPoint3 p; e->pnt(u, v, w, p);
inv3x3(jac, invjac);
e->getGradShapeFunctions(u, v, w, grads);
- for (int j = 0; j < nbNodes; j++){
+ for(int j = 0; j < nbSF; j++){
Grads[j] = SVector3(invjac[0][0] * grads[j][0] + invjac[0][1] * grads[j][1] +
invjac[0][2] * grads[j][2],
invjac[1][0] * grads[j][0] + invjac[1][1] * grads[j][1] +
@@ -58,26 +57,26 @@ class orthogonalTerm : public helmholtzTerm<double> {
invjac[2][2] * grads[j][2]);
}
SVector3 N (jac[2][0], jac[2][1], jac[2][2]);
- for (int j = 0; j < nbNodes; j++)
- for (int k = 0; k <= j; k++)
+ for(int j = 0; j < nbSF; j++)
+ for(int k = 0; k <= j; k++)
mat(j, k) += dot(crossprod(Grads[j], Grads[k]), N) * weight * detJ;
}
- for (int j = 0; j < nbNodes; j++)
- for (int k = 0; k < j; k++)
+ for(int j = 0; j < nbSF; j++)
+ for(int k = 0; k < j; k++)
mat(k, j) = -1.* mat(j, k);
//2) compute vector m(i) = mat(i,j)*val(j)
- fullVector<double> val(nbNodes);
+ fullVector<double> val(nbSF);
val.scale(0.);
- for (int i = 0; i < nbNodes; i++){
- std::map<MVertex*, double>::iterator it = _distance_map->find(e->getVertex(i));
+ for(int i = 0; i < nbSF; i++){
+ std::map<MVertex*, double>::iterator it = _distance_map->find(e->getShapeFunctionNode(i));
val(i) = it->second;
}
/* if( withDof){ */
- /* for (int i = 0; i < nbNodes; i++){ */
- /* _dofView->getDofValue( e->getVertex(i), 0, 1, val(i)); */
+ /* for (int i = 0; i < nbSF; i++){ */
+ /* _dofView->getDofValue( e->getShapeFunctionNode(i), 0, 1, val(i)); */
/* printf("val=%g \n", val(i)); */
/* } */
/* } */
@@ -85,18 +84,16 @@ class orthogonalTerm : public helmholtzTerm<double> {
/* printf("with no dof \n"); */
/* exit(1); */
/* /\* PViewData *data = view->getData(); *\/ */
- /* /\* for (int i = 0; i < nbNodes; i++){ *\/ */
- /* /\* data->getValue(0, e->, e->getNum(), e->getVertex(i)->getNum(), nod, 0, val(i)); *\/ */
+ /* /\* for (int i = 0; i < nbSF; i++){ *\/ */
+ /* /\* data->getValue(0, e->, e->getNum(), e->getShapeFunctionNode(i)->getNum(), nod, 0, val(i)); *\/ */
/* /\* } *\/ */
/* } */
m.scale(0.);
- for (int i = 0; i < nbNodes; i++)
- for (int j = 0; j < nbNodes; j++)
- m(i) += -mat(i,j)*val(j);
-
+ for(int i = 0; i < nbSF; i++)
+ for(int j = 0; j < nbSF; j++)
+ m(i) += -mat(i,j) * val(j);
}
-
};
#endif
diff --git a/Solver/terms.h b/Solver/terms.h
index 7d52068f31..a336c2490a 100644
--- a/Solver/terms.h
+++ b/Solver/terms.h
@@ -26,7 +26,7 @@ class BilinearTermBase
{
public :
virtual ~BilinearTermBase() {}
- virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m) = 0 ;
+ virtual void get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) = 0 ;
};
template<class T1,class T2> class BilinearTerm : public BilinearTermBase
@@ -35,7 +35,7 @@ template<class T1,class T2> class BilinearTerm : public BilinearTermBase
FunctionSpace<T1>& space1;
FunctionSpace<T2>& space2;
public :
- BilinearTerm(FunctionSpace<T1>& space1_,FunctionSpace<T2>& space2_) : space1(space1_),space2(space2_) {}
+ BilinearTerm(FunctionSpace<T1>& space1_, FunctionSpace<T2>& space2_) : space1(space1_), space2(space2_) {}
virtual ~BilinearTerm() {}
};
@@ -43,7 +43,7 @@ class LinearTermBase
{
public:
virtual ~LinearTermBase() {}
- virtual void get(MElement *ele,int npts,IntPt *GP,fullVector<double> &v) =0;
+ virtual void get(MElement *ele, int npts, IntPt *GP, fullVector<double> &v) = 0;
};
template<class T1> class LinearTerm : public LinearTermBase
@@ -59,7 +59,7 @@ class ScalarTermBase
{
public :
virtual ~ScalarTermBase() {}
- virtual void get(MElement *ele,int npts,IntPt *GP,double &val) =0;
+ virtual void get(MElement *ele, int npts, IntPt *GP, double &val) = 0;
};
class ScalarTerm : public ScalarTermBase
@@ -68,104 +68,91 @@ class ScalarTerm : public ScalarTermBase
virtual ~ScalarTerm() {}
};
-
-template<class T1,class T2> class BilinearTermToScalarTerm : public ScalarTerm
+template<class T1, class T2> class BilinearTermToScalarTerm : public ScalarTerm
{
- BilinearTerm<T1,T2> &bilterm;
+ BilinearTerm<T1, T2> &bilterm;
public :
- BilinearTermToScalarTerm(BilinearTerm<T1,T2> &bilterm_): bilterm(bilterm_){}
+ BilinearTermToScalarTerm(BilinearTerm<T1, T2> &bilterm_): bilterm(bilterm_){}
virtual ~BilinearTermToScalarTerm() {}
- virtual void get(MElement *ele,int npts,IntPt *GP,double &val)
+ virtual void get(MElement *ele, int npts, IntPt *GP, double &val)
{
fullMatrix<double> localMatrix;
- bilterm.get(ele,npts,GP,localMatrix);
- val=localMatrix(0,0);
+ bilterm.get(ele, npts, GP, localMatrix);
+ val = localMatrix(0, 0);
}
};
-
class ScalarTermConstant : public ScalarTerm
{
double val;
public :
- ScalarTermConstant(double val_=1.0): val(val_) {}
+ ScalarTermConstant(double val_ = 1.0) : val(val_) {}
virtual ~ScalarTermConstant() {}
- virtual void get(MElement *ele,int npts,IntPt *GP,double &val)
+ virtual void get(MElement *ele, int npts, IntPt *GP, double &val)
{
double jac[3][3];
- val=0;
- for (int i = 0; i < npts; i++)
- {
- const double u = GP[i].pt[0];const double v = GP[i].pt[1];const double w = GP[i].pt[2];
- const double weight = GP[i].weight;const double detJ = ele->getJacobian(u, v, w, jac);
- val+=weight*detJ;
+ val = 0;
+ for(int i = 0; i < npts; i++){
+ const double u = GP[i].pt[0]; const double v = GP[i].pt[1]; const double w = GP[i].pt[2];
+ const double weight = GP[i].weight; const double detJ = ele->getJacobian(u, v, w, jac);
+ val += weight * detJ;
}
}
- virtual void get(MVertex *ver,double &val)
+ virtual void get(MVertex *ver, double &val)
{
- val=1;
+ val = 1;
}
};
-
-
-template<class T1,class T2> class LaplaceTerm : public BilinearTerm<T1,T2>
+template<class T1, class T2> class LaplaceTerm : public BilinearTerm<T1, T2>
{
public :
- LaplaceTerm(FunctionSpace<T1>& space1_,FunctionSpace<T2>& space2_) : BilinearTerm<T1,T2>(space1_,space2_)
+ LaplaceTerm(FunctionSpace<T1>& space1_, FunctionSpace<T2>& space2_) : BilinearTerm<T1, T2>(space1_, space2_)
{}
virtual ~LaplaceTerm() {}
- virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m)
+ virtual void get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m)
{
- Msg::Error("LaplaceTerm<S1,S2> w/ S1!=S2 not implemented");
+ Msg::Error("LaplaceTerm<S1, S2> w/ S1 != S2 not implemented");
}
- virtual void get(MVertex *ver,fullMatrix<double> &m)
+ virtual void get(MVertex *ver, fullMatrix<double> &m)
{
- Msg::Error("LaplaceTerm<S1,S2> w/ S1!=S2 not implemented");
+ Msg::Error("LaplaceTerm<S1, S2> w/ S1 != S2 not implemented");
}
}; // class
-
-
-template<class T1> class LaplaceTerm<T1,T1> : public BilinearTerm<T1,T1> // symmetric
+template<class T1> class LaplaceTerm<T1, T1> : public BilinearTerm<T1, T1> // symmetric
{
double diffusivity;
public :
- LaplaceTerm(FunctionSpace<T1>& space1_, double diff=1) : BilinearTerm<T1,T1>(space1_,space1_), diffusivity(diff)
- {}
+ LaplaceTerm(FunctionSpace<T1>& space1_, double diff = 1) :
+ BilinearTerm<T1, T1>(space1_, space1_), diffusivity(diff) {}
virtual ~LaplaceTerm() {}
- virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m)
+ virtual void get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m)
{
- int nbFF = BilinearTerm<T1,T1>::space1.getNumKeys(ele);
+ int nbFF = BilinearTerm<T1, T1>::space1.getNumKeys(ele);
double jac[3][3];
m.resize(nbFF, nbFF);
m.setAll(0.);
- for (int i = 0; i < npts; i++)
- {
+ for(int i = 0; i < npts; i++){
const double u = GP[i].pt[0]; const double v = GP[i].pt[1]; const double w = GP[i].pt[2];
const double weight = GP[i].weight; const double detJ = ele->getJacobian(u, v, w, jac);
std::vector<typename TensorialTraits<T1>::GradType> Grads;
- BilinearTerm<T1,T1>::space1.gradf(ele,u, v, w, Grads);
- for (int j = 0; j < nbFF; j++)
- {
- for (int k = j; k < nbFF; k++)
- {
- double contrib=weight * detJ * dot(Grads[j],Grads[k]) * diffusivity;
- m(j,k)+=contrib;
- if (j!=k) m(k,j)+=contrib;
+ BilinearTerm<T1, T1>::space1.gradf(ele, u, v, w, Grads);
+ for(int j = 0; j < nbFF; j++){
+ for(int k = j; k < nbFF; k++){
+ double contrib = weight * detJ * dot(Grads[j], Grads[k]) * diffusivity;
+ m(j, k) += contrib;
+ if(j != k) m(k, j) += contrib;
}
}
}
-// m.print("");
-// exit(0);
}
}; // class
-
class IsotropicElasticTerm : public BilinearTerm<SVector3,SVector3>
{
protected :
- double E,nu;
+ double E, nu;
bool sym;
fullMatrix<double> H;/* =
{ {C11, C12, C12, 0, 0, 0},
@@ -174,75 +161,69 @@ class IsotropicElasticTerm : public BilinearTerm<SVector3,SVector3>
{ 0, 0, 0, C44, 0, 0},
{ 0, 0, 0, 0, C44, 0},
{ 0, 0, 0, 0, 0, C44} };*/
-
public :
- IsotropicElasticTerm(FunctionSpace<SVector3>& space1_,FunctionSpace<SVector3>& space2_,double E_,double nu_) : BilinearTerm<SVector3,SVector3>(space1_,space2_),E(E_),nu(nu_),H(6,6)
+ IsotropicElasticTerm(FunctionSpace<SVector3>& space1_, FunctionSpace<SVector3>& space2_, double E_, double nu_) :
+ BilinearTerm<SVector3, SVector3>(space1_, space2_), E(E_), nu(nu_), H(6, 6)
{
double FACT = E / (1 + nu);
double C11 = FACT * (1 - nu) / (1 - 2 * nu);
double C12 = FACT * nu / (1 - 2 * nu);
double C44 = (C11 - C12) / 2;
H.scale(0.);
- for (int i=0;i<3;++i) {H(i,i)=C11;H(i+3,i+3)=C44;}
- H(1,0)=H(0,1)=H(2,0)=H(0,2)=H(1,2)=H(2,1)=C12;
- sym=(&space1_==&space2_);
+ for(int i = 0; i < 3; ++i) { H(i, i) = C11; H(i + 3, i + 3) = C44; }
+ H(1, 0) = H(0, 1) = H(2, 0) = H(0, 2) = H(1, 2) = H(2, 1) = C12;
+ sym = (&space1_ == &space2_);
}
- IsotropicElasticTerm(FunctionSpace<SVector3>& space1_,double E_,double nu_) : BilinearTerm<SVector3,SVector3>(space1_,space1_),E(E_),nu(nu_),H(6,6)
+ IsotropicElasticTerm(FunctionSpace<SVector3>& space1_, double E_, double nu_) :
+ BilinearTerm<SVector3, SVector3>(space1_, space1_), E(E_), nu(nu_), H(6, 6)
{
double FACT = E / (1 + nu);
double C11 = FACT * (1 - nu) / (1 - 2 * nu);
double C12 = FACT * nu / (1 - 2 * nu);
double C44 = (C11 - C12) / 2;
- FACT = E / (1-nu*nu);
- C11 = FACT;
+ FACT = E / (1 - nu * nu);
+ C11 = FACT;
C12 = nu * FACT;
- C44 = (1.-nu)*.5*FACT;
-
+ C44 = (1. - nu) * .5 * FACT;
H.scale(0.);
- for (int i=0;i<3;++i) {H(i,i)=C11;H(i+3,i+3)=C44;}
- H(1,0)=H(0,1)=H(2,0)=H(0,2)=H(1,2)=H(2,1)=C12;
- sym=true;
+ for(int i = 0; i < 3; ++i) { H(i, i) = C11; H(i + 3, i + 3) = C44; }
+ H(1, 0) = H(0, 1) = H(2, 0) = H(0, 2) = H(1, 2) = H(2, 1) = C12;
+ sym = true;
}
virtual ~IsotropicElasticTerm() {}
- virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m)
+ virtual void get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m)
{
- if (ele->getParent()) ele=ele->getParent();
- if (sym)
- {
- int nbFF = BilinearTerm<SVector3,SVector3>::space1.getNumKeys(ele);
+ if(ele->getParent()) ele = ele->getParent();
+ if(sym){
+ int nbFF = BilinearTerm<SVector3, SVector3>::space1.getNumKeys(ele);
double jac[3][3];
fullMatrix<double> B(6, nbFF);
fullMatrix<double> BTH(nbFF, 6);
fullMatrix<double> BT(nbFF, 6);
- //printf("npts : %d\n",npts);
m.resize(nbFF, nbFF);
m.setAll(0.);
//std::cout << m.size1() << " " << m.size2() << std::endl;
- for (int i = 0; i < npts; i++)
- {
+ for(int i = 0; i < npts; i++){
const double u = GP[i].pt[0]; const double v = GP[i].pt[1]; const double w = GP[i].pt[2];
const double weight = GP[i].weight; const double detJ = ele->getJacobian(u, v, w, jac);
std::vector<TensorialTraits<SVector3>::GradType> Grads;
- BilinearTerm<SVector3,SVector3>::space1.gradf(ele,u, v, w, Grads); // a optimiser ??
- for (int j = 0; j < nbFF; j++)
- {
- BT(j, 0) = B(0, j) = Grads[j](0,0);
- BT(j, 1) = B(1, j) = Grads[j](1,1);
- BT(j, 2) = B(2, j) = Grads[j](2,2);
- BT(j, 3) = B(3, j) = Grads[j](0,1)+Grads[j](1,0);
- BT(j, 4) = B(4, j) = Grads[j](1,2)+Grads[j](2,1);
- BT(j, 5) = B(5, j) = Grads[j](0,2)+Grads[j](2,0);
+ BilinearTerm<SVector3, SVector3>::space1.gradf(ele, u, v, w, Grads); // a optimiser ??
+ for(int j = 0; j < nbFF; j++){
+ BT(j, 0) = B(0, j) = Grads[j](0, 0);
+ BT(j, 1) = B(1, j) = Grads[j](1, 1);
+ BT(j, 2) = B(2, j) = Grads[j](2, 2);
+ BT(j, 3) = B(3, j) = Grads[j](0, 1) + Grads[j](1, 0);
+ BT(j, 4) = B(4, j) = Grads[j](1, 2) + Grads[j](2, 1);
+ BT(j, 5) = B(5, j) = Grads[j](0, 2) + Grads[j](2, 0);
}
BTH.setAll(0.);
-
BTH.gemm(BT, H);
m.gemm(BTH, B, weight * detJ, 1.);
}
}
- else
- {
- int nbFF1 = BilinearTerm<SVector3,SVector3>::space1.getNumKeys(ele);
- int nbFF2 = BilinearTerm<SVector3,SVector3>::space2.getNumKeys(ele);
+ else{
+ int nbFF1 = BilinearTerm<SVector3, SVector3>::space1.getNumKeys(ele);
+ int nbFF2 = BilinearTerm<SVector3, SVector3>::space2.getNumKeys(ele);
double jac[3][3];
fullMatrix<double> B(6, nbFF2);
fullMatrix<double> BTH(nbFF2, 6);
@@ -250,31 +231,28 @@ class IsotropicElasticTerm : public BilinearTerm<SVector3,SVector3>
m.resize(nbFF1, nbFF2);
m.setAll(0.);
// Sum on Gauss Points i
- for (int i = 0; i < npts; i++)
- {
+ for(int i = 0; i < npts; i++){
const double u = GP[i].pt[0]; const double v = GP[i].pt[1]; const double w = GP[i].pt[2];
const double weight = GP[i].weight; const double detJ = ele->getJacobian(u, v, w, jac);
std::vector<TensorialTraits<SVector3>::GradType> Grads;// tableau de matrices...
std::vector<TensorialTraits<SVector3>::GradType> GradsT;// tableau de matrices...
- BilinearTerm<SVector3,SVector3>::space1.gradf(ele,u, v, w, Grads);
- BilinearTerm<SVector3,SVector3>::space2.gradf(ele,u, v, w, GradsT);
- for (int j = 0; j < nbFF1; j++)
- {
- BT(j, 0) = Grads[j](0,0);
- BT(j, 1) = Grads[j](1,1);
- BT(j, 2) = Grads[j](2,2);
- BT(j, 3) = Grads[j](0,1)+Grads[j](1,0);
- BT(j, 4) = Grads[j](1,2)+Grads[j](2,1);
- BT(j, 5) = Grads[j](0,2)+Grads[j](2,0);
+ BilinearTerm<SVector3, SVector3>::space1.gradf(ele, u, v, w, Grads);
+ BilinearTerm<SVector3, SVector3>::space2.gradf(ele, u, v, w, GradsT);
+ for(int j = 0; j < nbFF1; j++){
+ BT(j, 0) = Grads[j](0, 0);
+ BT(j, 1) = Grads[j](1, 1);
+ BT(j, 2) = Grads[j](2, 2);
+ BT(j, 3) = Grads[j](0, 1) + Grads[j](1, 0);
+ BT(j, 4) = Grads[j](1, 2) + Grads[j](2, 1);
+ BT(j, 5) = Grads[j](0, 2) + Grads[j](2, 0);
}
- for (int j = 0; j < nbFF2; j++)
- {
- B(0, j) = GradsT[j](0,0);
- B(1, j) = GradsT[j](1,1);
- B(2, j) = GradsT[j](2,2);
- B(3, j) = GradsT[j](0,1)+GradsT[j](1,0);
- B(4, j) = GradsT[j](1,2)+GradsT[j](2,1);
- B(5, j) = GradsT[j](0,2)+GradsT[j](2,0);
+ for(int j = 0; j < nbFF2; j++){
+ B(0, j) = GradsT[j](0, 0);
+ B(1, j) = GradsT[j](1, 1);
+ B(2, j) = GradsT[j](2, 2);
+ B(3, j) = GradsT[j](0, 1) + GradsT[j](1, 0);
+ B(4, j) = GradsT[j](1, 2) + GradsT[j](2, 1);
+ B(5, j) = GradsT[j](0, 2) + GradsT[j](2, 0);
}
BTH.setAll(0.);
BTH.gemm(BT, H);
@@ -286,63 +264,66 @@ class IsotropicElasticTerm : public BilinearTerm<SVector3,SVector3>
}; // class
inline double dot(const double &a, const double &b)
-{ return a*b; }
+{
+ return a * b;
+}
template<class T1> class LoadTerm : public LinearTerm<T1>
{
simpleFunction<typename TensorialTraits<T1>::ValType> &Load;
public :
- LoadTerm(FunctionSpace<T1>& space1_,simpleFunction<typename TensorialTraits<T1>::ValType> &Load_) :LinearTerm<T1>(space1_),Load(Load_) {}
+ LoadTerm(FunctionSpace<T1>& space1_, simpleFunction<typename TensorialTraits<T1>::ValType> &Load_) :
+ LinearTerm<T1>(space1_), Load(Load_) {}
virtual ~LoadTerm() {}
- virtual void get(MElement *ele,int npts,IntPt *GP,fullVector<double> &m)
+ virtual void get(MElement *ele, int npts, IntPt *GP, fullVector<double> &m)
{
- if (ele->getParent()) ele=ele->getParent();
- int nbFF=LinearTerm<T1>::space1.getNumKeys(ele);
+ if(ele->getParent()) ele = ele->getParent();
+ int nbFF = LinearTerm<T1>::space1.getNumKeys(ele);
double jac[3][3];
m.resize(nbFF);
m.scale(0.);
- for (int i = 0; i < npts; i++)
- {
+ for(int i = 0; i < npts; i++){
const double u = GP[i].pt[0]; const double v = GP[i].pt[1]; const double w = GP[i].pt[2];
const double weight = GP[i].weight; const double detJ = ele->getJacobian(u, v, w, jac);
std::vector<typename TensorialTraits<T1>::ValType> Vals;
LinearTerm<T1>::space1.f(ele, u, v, w, Vals);
SPoint3 p;
ele->pnt(u, v, w, p);
- typename TensorialTraits<T1>::ValType load=Load(p.x(),p.y(),p.z());
- for (int j = 0; j < nbFF ; ++j)
- {
+ typename TensorialTraits<T1>::ValType load = Load(p.x(), p.y(), p.z());
+ for(int j = 0; j < nbFF ; ++j){
m(j) += dot(Vals[j], load) * weight * detJ;
}
}
}
};
-
-class LagrangeMultiplierTerm : public BilinearTerm<SVector3,double>
+class LagrangeMultiplierTerm : public BilinearTerm<SVector3, double>
{
SVector3 _d;
public :
LagrangeMultiplierTerm(FunctionSpace<SVector3>& space1_, FunctionSpace<double>& space2_, const SVector3 &d) :
- BilinearTerm<SVector3,double>(space1_, space2_) {for(int i=0; i < 3; i++) _d(i) = d(i);}
+ BilinearTerm<SVector3, double>(space1_, space2_)
+ {
+ for(int i = 0; i < 3; i++) _d(i) = d(i);
+ }
virtual ~LagrangeMultiplierTerm() {}
virtual void get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m)
{
- int nbFF1 = BilinearTerm<SVector3,double>::space1.getNumKeys(ele); //nbVertices*nbcomp of parent
- int nbFF2 = BilinearTerm<SVector3,double>::space2.getNumKeys(ele); //nbVertices of boundary
+ int nbFF1 = BilinearTerm<SVector3, double>::space1.getNumKeys(ele); //nbVertices*nbcomp of parent
+ int nbFF2 = BilinearTerm<SVector3, double>::space2.getNumKeys(ele); //nbVertices of boundary
double jac[3][3];
m.resize(nbFF1, nbFF2);
m.setAll(0.);
- for (int i = 0; i < npts; i++) {
+ for(int i = 0; i < npts; i++){
double u = GP[i].pt[0]; double v = GP[i].pt[1]; double w = GP[i].pt[2];
const double weight = GP[i].weight; const double detJ = ele->getJacobian(u, v, w, jac);
std::vector<TensorialTraits<SVector3>::ValType> Vals;
std::vector<TensorialTraits<double>::ValType> ValsT;
BilinearTerm<SVector3,double>::space1.f(ele, u, v, w, Vals);
BilinearTerm<SVector3,double>::space2.f(ele, u, v, w, ValsT);
- for (int j = 0; j < nbFF1; j++) {
- for (int k = 0; k < nbFF2; k++) {
+ for(int j = 0; j < nbFF1; j++){
+ for(int k = 0; k < nbFF2; k++){
m(j, k) += dot(Vals[j], _d) * ValsT[k] * weight * detJ;
}
}
@@ -350,79 +331,67 @@ class LagrangeMultiplierTerm : public BilinearTerm<SVector3,double>
}
};
-
class LagMultTerm : public BilinearTerm<SVector3, SVector3>
{
-
private :
-
double _eqfac;
-
public :
LagMultTerm(FunctionSpace<SVector3>& space1_, FunctionSpace<SVector3>& space2_, double eqfac = 1.0) :
- BilinearTerm<SVector3,SVector3>(space1_, space2_), _eqfac(eqfac) {;}
+ BilinearTerm<SVector3, SVector3>(space1_, space2_), _eqfac(eqfac) {}
virtual ~LagMultTerm() {}
virtual void get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m)
{
- int nbFF1 = BilinearTerm<SVector3,SVector3>::space1.getNumKeys(ele); //nbVertices*nbcomp of parent
- int nbFF2 = BilinearTerm<SVector3,SVector3>::space2.getNumKeys(ele); //nbVertices of boundary
+ int nbFF1 = BilinearTerm<SVector3, SVector3>::space1.getNumKeys(ele); //nbVertices*nbcomp of parent
+ int nbFF2 = BilinearTerm<SVector3, SVector3>::space2.getNumKeys(ele); //nbVertices of boundary
double jac[3][3];
m.resize(nbFF1, nbFF2);
m.setAll(0.);
- for (int i = 0; i < npts; i++) {
+ for(int i = 0; i < npts; i++) {
double u = GP[i].pt[0]; double v = GP[i].pt[1]; double w = GP[i].pt[2];
const double weight = GP[i].weight; const double detJ = ele->getJacobian(u, v, w, jac);
std::vector<TensorialTraits<SVector3>::ValType> Vals;
std::vector<TensorialTraits<SVector3>::ValType> ValsT;
BilinearTerm<SVector3,SVector3>::space1.f(ele, u, v, w, Vals);
BilinearTerm<SVector3,SVector3>::space2.f(ele, u, v, w, ValsT);
- for (int j = 0; j < nbFF1; j++) {
- for (int k = 0; k < nbFF2; k++) {
+ for(int j = 0; j < nbFF1; j++){
+ for(int k = 0; k < nbFF2; k++){
m(j, k) += _eqfac * dot(Vals[j], ValsT[k]) * weight * detJ;
}
}
}
- //m.print();
}
};
-
template<class T1> class LoadTermOnBorder : public LinearTerm<T1>
{
private :
-
double _eqfac;
simpleFunction<typename TensorialTraits<T1>::ValType> &Load;
-
public :
-
- LoadTermOnBorder(FunctionSpace<T1>& space1_,simpleFunction<typename TensorialTraits<T1>::ValType> &Load_, double eqfac = 1.0) : LinearTerm<T1>(space1_), _eqfac(eqfac), Load(Load_) {}
+ LoadTermOnBorder(FunctionSpace<T1>& space1_, simpleFunction<typename TensorialTraits<T1>::ValType> &Load_, double eqfac = 1.0) :
+ LinearTerm<T1>(space1_), _eqfac(eqfac), Load(Load_) {}
virtual ~LoadTermOnBorder() {}
-
- virtual void get(MElement *ele,int npts,IntPt *GP,fullVector<double> &m)
+ virtual void get(MElement *ele, int npts, IntPt *GP, fullVector<double> &m)
{
- MElement * elep;
- if (ele->getParent()) elep=ele->getParent();
- int nbFF=LinearTerm<T1>::space1.getNumKeys(ele);
+ MElement *elep;
+ if (ele->getParent()) elep = ele->getParent();
+ int nbFF = LinearTerm<T1>::space1.getNumKeys(ele);
double jac[3][3];
m.resize(nbFF);
m.scale(0.);
- for (int i = 0; i < npts; i++)
- {
+ for(int i = 0; i < npts; i++){
const double u = GP[i].pt[0]; const double v = GP[i].pt[1]; const double w = GP[i].pt[2];
const double weight = GP[i].weight; const double detJ = ele->getJacobian(u, v, w, jac);
std::vector<typename TensorialTraits<T1>::ValType> Vals;
LinearTerm<T1>::space1.f(ele, u, v, w, Vals);
SPoint3 p;
ele->pnt(u, v, w, p);
- typename TensorialTraits<T1>::ValType load=Load(p.x(),p.y(),p.z());
- for (int j = 0; j < nbFF ; ++j)
- {
- m(j) += _eqfac*dot(Vals[j], load) * weight * detJ;
+ typename TensorialTraits<T1>::ValType load = Load(p.x(), p.y(), p.z());
+ for(int j = 0; j < nbFF ; ++j){
+ m(j) += _eqfac * dot(Vals[j], load) * weight * detJ;
}
}
}
};
-
#endif// _TERMS_H_
--
GitLab