diff --git a/Geo/GEntity.h b/Geo/GEntity.h
index 9887aa4a7a149cd2e7052f99163cd98299db16fe..668ee76d4d7d7feb0240370443864c0d84a5fe0b 100644
--- a/Geo/GEntity.h
+++ b/Geo/GEntity.h
@@ -222,14 +222,15 @@ class GEntity {
void setTag(int tag) { _tag = tag; }
// get/set physical entities
- virtual void addPhysicalEntity(int physicalTag) {
+ virtual void addPhysicalEntity(int physicalTag)
+ {
physicals.push_back(physicalTag);
}
- virtual std::vector<int> getPhysicalEntities() {
+ virtual std::vector<int> getPhysicalEntities()
+ {
return physicals;
}
-
// returns the tag of the entity that its master entity (for mesh)
int meshMaster() const { return _meshMaster; }
void setMeshMaster(int m) { _meshMaster = m; }
diff --git a/Geo/GModelFactory.cpp b/Geo/GModelFactory.cpp
index a6bfb66a5511fca5abb1a1d1b3ee92eb56c58663..216f9dbe273c6384cf527456280875b4d94daffd 100644
--- a/Geo/GModelFactory.cpp
+++ b/Geo/GModelFactory.cpp
@@ -15,10 +15,8 @@
#include "MLine.h"
#include "GModel.h"
-
GVertex *GeoFactory::addVertex(GModel *gm, double x, double y, double z, double lc)
{
-
int num = gm->getMaxElementaryNumber(0) + 1;
x *= CTX::instance()->geom.scalingFactor;
@@ -36,12 +34,10 @@ GVertex *GeoFactory::addVertex(GModel *gm, double x, double y, double z, double
gm->add(v);
return v;
-
}
GEdge *GeoFactory::addLine(GModel *gm, GVertex *start, GVertex *end)
{
-
int num = gm->getMaxElementaryNumber(1) + 1;
List_T *iList = List_Create(2, 2, sizeof(int));
int tagBeg = start->tag();
@@ -61,33 +57,32 @@ GEdge *GeoFactory::addLine(GModel *gm, GVertex *start, GVertex *end)
gm->add(e);
return e;
-
}
GFace *GeoFactory::addPlanarFace(GModel *gm, std::vector< std::vector<GEdge *> > edges)
{
//create line loops
- std::vector<EdgeLoop *> vecLoops;
- int nLoops = edges.size();
- for (int i=0; i< nLoops; i++){
- int numl = gm->getMaxElementaryNumber(1) + i;
- while (FindEdgeLoop(numl)){
- numl++;
- if (!FindEdgeLoop(numl)) break;
- }
- int nl=(int)edges[i].size();
- List_T *iListl = List_Create(nl, nl, sizeof(int));
- for(int j = 0; j < nl; j++){
- int numEdge = edges[i][j]->tag();
- List_Add(iListl, &numEdge);
- }
- sortEdgesInLoop(numl, iListl);
- EdgeLoop *l = Create_EdgeLoop(numl, iListl);
- vecLoops.push_back(l);
- Tree_Add(GModel::current()->getGEOInternals()->EdgeLoops, &l);
- l->Num = numl;
- List_Delete(iListl);
- }
+ std::vector<EdgeLoop *> vecLoops;
+ int nLoops = edges.size();
+ for (int i=0; i< nLoops; i++){
+ int numl = gm->getMaxElementaryNumber(1) + i;
+ while (FindEdgeLoop(numl)){
+ numl++;
+ if (!FindEdgeLoop(numl)) break;
+ }
+ int nl=(int)edges[i].size();
+ List_T *iListl = List_Create(nl, nl, sizeof(int));
+ for(int j = 0; j < nl; j++){
+ int numEdge = edges[i][j]->tag();
+ List_Add(iListl, &numEdge);
+ }
+ sortEdgesInLoop(numl, iListl);
+ EdgeLoop *l = Create_EdgeLoop(numl, iListl);
+ vecLoops.push_back(l);
+ Tree_Add(GModel::current()->getGEOInternals()->EdgeLoops, &l);
+ l->Num = numl;
+ List_Delete(iListl);
+ }
//create plane surfaces
int numf = gm->getMaxElementaryNumber(2) + 1;
@@ -109,55 +104,52 @@ GFace *GeoFactory::addPlanarFace(GModel *gm, std::vector< std::vector<GEdge *> >
gm->add(gf);
return gf;
-
}
GRegion* GeoFactory::addVolume (GModel *gm, std::vector<std::vector<GFace *> > faces)
{
-
//surface loop
- std::vector<SurfaceLoop *> vecLoops;
- int nLoops = faces.size();
- for (int i=0; i< nLoops; i++){
- int numfl = gm->getMaxElementaryNumber(2) + 1;
- while (FindSurfaceLoop(numfl)){
- numfl++;
- if (!FindSurfaceLoop(numfl)) break;
- }
- int nl=(int)faces[i].size();
- List_T *iListl = List_Create(nl, nl, sizeof(int));
- for(int j = 0; j < nl; j++){
- int numFace = faces[i][j]->tag();
- List_Add(iListl, &numFace);
- }
- SurfaceLoop *l = Create_SurfaceLoop(numfl, iListl);
- vecLoops.push_back(l);
- Tree_Add(GModel::current()->getGEOInternals()->SurfaceLoops, &l);
- List_Delete(iListl);
- }
-
- //volume
- int numv = gm->getMaxElementaryNumber(3) + 1;
- Volume *v = Create_Volume(numv, MSH_VOLUME);
- List_T *iList = List_Create(nLoops, nLoops, sizeof(int));
- for (int i=0; i< vecLoops.size(); i++){
- int numl = vecLoops[i]->Num;
- List_Add(iList, &numl);
- }
- setVolumeSurfaces(v, iList);
- List_Delete(iList);
- Tree_Add(GModel::current()->getGEOInternals()->Volumes, &v);
- v->Typ = MSH_VOLUME;
- v->Num = numv;
-
- //gmsh volume
+ std::vector<SurfaceLoop *> vecLoops;
+ int nLoops = faces.size();
+ for (int i=0; i< nLoops; i++){
+ int numfl = gm->getMaxElementaryNumber(2) + 1;
+ while (FindSurfaceLoop(numfl)){
+ numfl++;
+ if (!FindSurfaceLoop(numfl)) break;
+ }
+ int nl=(int)faces[i].size();
+ List_T *iListl = List_Create(nl, nl, sizeof(int));
+ for(int j = 0; j < nl; j++){
+ int numFace = faces[i][j]->tag();
+ List_Add(iListl, &numFace);
+ }
+ SurfaceLoop *l = Create_SurfaceLoop(numfl, iListl);
+ vecLoops.push_back(l);
+ Tree_Add(GModel::current()->getGEOInternals()->SurfaceLoops, &l);
+ List_Delete(iListl);
+ }
+
+ //volume
+ int numv = gm->getMaxElementaryNumber(3) + 1;
+ Volume *v = Create_Volume(numv, MSH_VOLUME);
+ List_T *iList = List_Create(nLoops, nLoops, sizeof(int));
+ for (int i=0; i< vecLoops.size(); i++){
+ int numl = vecLoops[i]->Num;
+ List_Add(iList, &numl);
+ }
+ setVolumeSurfaces(v, iList);
+ List_Delete(iList);
+ Tree_Add(GModel::current()->getGEOInternals()->Volumes, &v);
+ v->Typ = MSH_VOLUME;
+ v->Num = numv;
+
+ //gmsh volume
GRegion *gr = new gmshRegion(gm,v);
gm->add(gr);
-
+
return gr;
}
-
//---------------------------------------------------------------------------------
#if defined(HAVE_OCC)
#include "OCCIncludes.h"
diff --git a/Geo/GRegion.cpp b/Geo/GRegion.cpp
index fa43377ae45a9a9e73005682f3c9b30954bdd9f2..bbf51bd3c9c944bbeee12278e2d9450a00cbd1a2 100644
--- a/Geo/GRegion.cpp
+++ b/Geo/GRegion.cpp
@@ -300,7 +300,8 @@ void GRegion::replaceFaces (std::list<GFace*> &new_faces)
}
double GRegion::computeSolidProperties (std::vector<double> cg,
- std::vector<double> inertia){
+ std::vector<double> inertia)
+{
std::list<GFace*>::iterator it = l_faces.begin();
std::list<int>::iterator itdir = l_dirs.begin();
double volumex = 0;
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 2f5dbc375e34c8288b4f1ea855a7aa8888b8c719..47e4c16e8b611ed00cc9efb1fd58ff19afbb748d 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -223,15 +223,15 @@ class MElement
int order=-1);
virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3],
int order=-1);
- //const fullMatrix<double> &getShapeFunctionsAtIntegrationPoints (int integrationOrder, int functionSpaceOrder=-1);
- const fullMatrix<double> &getGradShapeFunctionsAtIntegrationPoints (int integrationOrder, int functionSpaceOrder=-1);
+ const fullMatrix<double> &getGradShapeFunctionsAtIntegrationPoints
+ (int integrationOrder, int functionSpaceOrder=-1);
const fullMatrix<double> &getGradShapeFunctionsAtNodes (int functionSpaceOrder=-1);
+
// return the Jacobian of the element evaluated at point (u,v,w) in
// parametric coordinates
double getJacobian(const fullMatrix<double> &gsf, double jac[3][3]);
double getJacobian(double u, double v, double w, double jac[3][3]);
double getPrimaryJacobian(double u, double v, double w, double jac[3][3]);
- //bindings : double[3][3] is not easy to bind, we could use fullMatrix instead
double getJacobianDeterminant(double u, double v, double w);
// get the point in cartesian coordinates corresponding to the point
diff --git a/Geo/MPoint.h b/Geo/MPoint.h
index 993e7743358b532798fce65d8cb9a8395b2fe63d..92be1b850ba9a20a4c9cce729865a0a89bd80e3a 100644
--- a/Geo/MPoint.h
+++ b/Geo/MPoint.h
@@ -50,22 +50,27 @@ class MPoint : public MElement {
{
s[0][0] = s[0][1] = s[0][2] = 0.;
}
-
- virtual const polynomialBasis* getFunctionSpace(int o) const { return polynomialBases::find(MSH_PNT); }
- virtual const JacobianBasis* getJacobianFuncSpace(int o) const { return JacobianBases::find(MSH_PNT); }
+ virtual const polynomialBasis* getFunctionSpace(int o) const
+ {
+ return polynomialBases::find(MSH_PNT);
+ }
+ virtual const JacobianBasis* getJacobianFuncSpace(int o) const
+ {
+ return JacobianBases::find(MSH_PNT);
+ }
virtual bool isInside(double u, double v, double w)
{
return true;
}
virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
{
- static IntPt GQL[1];
- GQL[0].pt[0] = 0;
- GQL[0].pt[1] = 0;
- GQL[0].pt[2] = 0;
- GQL[0].weight = 1;
- *npts = 1;
- *pts = GQL;
+ static IntPt GQL[1];
+ GQL[0].pt[0] = 0;
+ GQL[0].pt[1] = 0;
+ GQL[0].pt[2] = 0;
+ GQL[0].weight = 1;
+ *npts = 1;
+ *pts = GQL;
}
};
diff --git a/Geo/MPrism.cpp b/Geo/MPrism.cpp
index 8333b2e576a4da9f860f6838d74a51012b573b86..d29fbe12dd88b38ab8a16d5c1d6071968ca88be3 100644
--- a/Geo/MPrism.cpp
+++ b/Geo/MPrism.cpp
@@ -119,28 +119,36 @@ void MPrism::getFaceInfo(const MFace &face, int &ithFace, int &sign, int &rot) c
}
else {
MVertex *v3 = _v[faces_prism(ithFace, 3)];
- if (v0 == face.getVertex(0) && v1 == face.getVertex(1) && v2 == face.getVertex(2) && v3 == face.getVertex(3)){
+ if (v0 == face.getVertex(0) && v1 == face.getVertex(1) &&
+ v2 == face.getVertex(2) && v3 == face.getVertex(3)){
sign = 1; rot = 0; return;
}
- if (v0 == face.getVertex(1) && v1 == face.getVertex(2) && v2 == face.getVertex(3) && v3 == face.getVertex(0)){
+ if (v0 == face.getVertex(1) && v1 == face.getVertex(2) &&
+ v2 == face.getVertex(3) && v3 == face.getVertex(0)){
sign = 1; rot = 1; return;
}
- if (v0 == face.getVertex(2) && v1 == face.getVertex(3) && v2 == face.getVertex(0) && v3 == face.getVertex(1)){
+ if (v0 == face.getVertex(2) && v1 == face.getVertex(3) &&
+ v2 == face.getVertex(0) && v3 == face.getVertex(1)){
sign = 1; rot = 2; return;
}
- if (v0 == face.getVertex(3) && v1 == face.getVertex(0) && v2 == face.getVertex(1) && v3 == face.getVertex(2)){
+ if (v0 == face.getVertex(3) && v1 == face.getVertex(0) &&
+ v2 == face.getVertex(1) && v3 == face.getVertex(2)){
sign = 1; rot = 3; return;
}
- if (v0 == face.getVertex(0) && v1 == face.getVertex(3) && v2 == face.getVertex(2) && v3 == face.getVertex(1)){
+ if (v0 == face.getVertex(0) && v1 == face.getVertex(3) &&
+ v2 == face.getVertex(2) && v3 == face.getVertex(1)){
sign = -1; rot = 0; return;
}
- if (v0 == face.getVertex(1) && v1 == face.getVertex(0) && v2 == face.getVertex(3) && v3 == face.getVertex(2)){
+ if (v0 == face.getVertex(1) && v1 == face.getVertex(0) &&
+ v2 == face.getVertex(3) && v3 == face.getVertex(2)){
sign = -1; rot = 1; return;
}
- if (v0 == face.getVertex(2) && v1 == face.getVertex(1) && v2 == face.getVertex(0) && v3 == face.getVertex(3)){
+ if (v0 == face.getVertex(2) && v1 == face.getVertex(1) &&
+ v2 == face.getVertex(0) && v3 == face.getVertex(3)){
sign = -1; rot = 2; return;
}
- if (v0 == face.getVertex(3) && v1 == face.getVertex(2) && v2 == face.getVertex(1) && v3 == face.getVertex(0)){
+ if (v0 == face.getVertex(3) && v1 == face.getVertex(2) &&
+ v2 == face.getVertex(1) && v3 == face.getVertex(0)){
sign = -1; rot = 3; return;
}
}
@@ -149,7 +157,8 @@ void MPrism::getFaceInfo(const MFace &face, int &ithFace, int &sign, int &rot) c
Msg::Error("Could not get face information for prism %d", getNum());
}
#include "Bindings.h"
-static MPrism18* MPrism18_binding(std::vector<MVertex*> v) {
+static MPrism18* MPrism18_binding(std::vector<MVertex*> v)
+{
return new MPrism18(v);
}
@@ -160,7 +169,8 @@ void MPrism::registerBindings(binding *b)
methodBinding *cm;
cm = cb->setConstructor<MPrism,MVertex*,MVertex*,MVertex*,MVertex*, MVertex*, MVertex*>();
cm->setArgNames("v0", "v1", "v2", "v3","v4","v5", NULL);
- cm->setDescription("Create a new prism with top triangle (v0,v1,v2) and bottom one (v3,v4,v5).");
+ cm->setDescription("Create a new prism with top triangle (v0,v1,v2) and "
+ "bottom one (v3,v4,v5).");
cm = cb->addMethod("getVolumeSign",&MPrism::getVolumeSign);
cm->setDescription("computes the sign of the element volume");
cm = cb->addMethod("revert",&MPrism::revert);
@@ -171,7 +181,7 @@ void MPrism::registerBindings(binding *b)
cb = b->addClass<MPrism18>("MPrism18");
cb->setDescription("A mesh second-order prism.");
cm = cb->addMethod("MPrism18",&MPrism18_binding);
-// cm = cb->setConstructor<MPrism18_binding,std::vector<MVertex*> >();
+ // cm = cb->setConstructor<MPrism18_binding,std::vector<MVertex*> >();
cm->setArgNames("vectorOfVertices", NULL);
cm->setDescription("Create a new prism with vertices in vectorV (length=18).");
cb->setParentClass<MPrism>();
diff --git a/Geo/MQuadrangle.h b/Geo/MQuadrangle.h
index d3c8e28b7b51c65f93126e39cca4442d9aa94a8b..2edbb12c7e07269c65d7c67555eef8b7e3a6b9fd 100644
--- a/Geo/MQuadrangle.h
+++ b/Geo/MQuadrangle.h
@@ -136,9 +136,10 @@ class MQuadrangle : public MElement {
virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
virtual double angleShapeMeasure();
virtual double distoShapeMeasure();
-// Computes the minimum inradius of the all the circles tangents to 3 of the 4
-// edges of the quad. If the 4 points of the quad are not planar, we compute
-// the mean plane due to the least-square criterion.
+ // Computes the minimum inradius of the all the circles tangents to
+ // 3 of the 4 edges of the quad. If the 4 points of the quad are not
+ // planar, we compute the mean plane due to the least-square
+ // criterion.
virtual double getInnerRadius();
private:
int edges_quad(const int edge, const int vert) const
@@ -377,18 +378,18 @@ class MQuadrangleN : public MQuadrangle {
}
virtual int getTypeForMSH() const
{
- if(_order==2 && _vs.size()+4==8) return MSH_QUA_8;
- if(_order==3 && _vs.size()+4==12) return MSH_QUA_12;
- if(_order==3 && _vs.size()+4==16) return MSH_QUA_16;
- if(_order==4 && _vs.size()+4==16) return MSH_QUA_16I;
- if(_order==4 && _vs.size()+4==25) return MSH_QUA_25;
- if(_order==5 && _vs.size()+4==20) return MSH_QUA_20;
- if(_order==5 && _vs.size()+4==36) return MSH_QUA_36;
- if(_order==6 && _vs.size()+4==49) return MSH_QUA_49;
- if(_order==7 && _vs.size()+4==64) return MSH_QUA_64;
- if(_order==8 && _vs.size()+4==81) return MSH_QUA_81;
- if(_order==9 && _vs.size()+4==100) return MSH_QUA_100;
- if(_order==10 && _vs.size()+4==121) return MSH_QUA_121;
+ if(_order==2 && _vs.size() + 4 == 8) return MSH_QUA_8;
+ if(_order==3 && _vs.size() + 4 == 12) return MSH_QUA_12;
+ if(_order==3 && _vs.size() + 4 == 16) return MSH_QUA_16;
+ if(_order==4 && _vs.size() + 4 == 16) return MSH_QUA_16I;
+ if(_order==4 && _vs.size() + 4 == 25) return MSH_QUA_25;
+ if(_order==5 && _vs.size() + 4 == 20) return MSH_QUA_20;
+ if(_order==5 && _vs.size() + 4 == 36) return MSH_QUA_36;
+ if(_order==6 && _vs.size() + 4 == 49) return MSH_QUA_49;
+ if(_order==7 && _vs.size() + 4 == 64) return MSH_QUA_64;
+ if(_order==8 && _vs.size() + 4 == 81) return MSH_QUA_81;
+ if(_order==9 && _vs.size() + 4 == 100) return MSH_QUA_100;
+ if(_order==10 && _vs.size() + 4 == 121) return MSH_QUA_121;
return 0;
}
virtual void revert()
@@ -402,30 +403,33 @@ class MQuadrangleN : public MQuadrangle {
};
template <class T>
-void inline sort2(T &a, T &b){
- if(b<a){
- T t=b;
- b=a;
- a=t;
+void inline sort2(T &a, T &b)
+{
+ if(b < a){
+ T t = b;
+ b = a;
+ a = t;
}
}
template <class T>
void sort4(T *t[3])
{
- sort2<T*>(t[0],t[1]);
- sort2<T*>(t[2],t[3]);
- sort2<T*>(t[0],t[2]);
- sort2<T*>(t[1],t[3]);
- sort2<T*>(t[1],t[2]);
+ sort2<T*>(t[0], t[1]);
+ sort2<T*>(t[2], t[3]);
+ sort2<T*>(t[0], t[2]);
+ sort2<T*>(t[1], t[3]);
+ sort2<T*>(t[1], t[2]);
}
struct compareMQuadrangleLexicographic
{
bool operator () (MQuadrangle *t1, MQuadrangle *t2) const
{
- MVertex *_v1[] = {t1->getVertex(0), t1->getVertex(1), t1->getVertex(2), t1->getVertex(3)};
- MVertex *_v2[] = {t2->getVertex(0), t2->getVertex(1), t2->getVertex(2), t2->getVertex(3)};
+ MVertex *_v1[] = {t1->getVertex(0), t1->getVertex(1),
+ t1->getVertex(2), t1->getVertex(3)};
+ MVertex *_v2[] = {t2->getVertex(0), t2->getVertex(1),
+ t2->getVertex(2), t2->getVertex(3)};
sort4(_v1);
sort4(_v2);
if(_v1[0] < _v2[0]) return true;
diff --git a/Geo/MTetrahedron.cpp b/Geo/MTetrahedron.cpp
index cc44a7c6fb460c2b9dca2c6bd4f01f5d1d0ccdc0..f93a07506650c2e85b7b334f09945993b2b5adca 100644
--- a/Geo/MTetrahedron.cpp
+++ b/Geo/MTetrahedron.cpp
@@ -364,9 +364,10 @@ void MTetrahedron::getFaceInfo(const MFace &face, int &ithFace, int &sign, int &
static std::vector<std::vector<int> > tetReverseIndices(20);
-const std::vector<int> &MTetrahedronN::_getReverseIndices (int order) {
- if(order>=tetReverseIndices.size())
- tetReverseIndices.resize(order+1);
+const std::vector<int> &MTetrahedronN::_getReverseIndices (int order)
+{
+ if(order >= tetReverseIndices.size())
+ tetReverseIndices.resize(order + 1);
std::vector<int> &r = tetReverseIndices[order];
if (r.size() != 0) return r;
//
diff --git a/Geo/MTetrahedron.h b/Geo/MTetrahedron.h
index fb735a6cc261e24c1dfd57e5d7cea837b13b64c6..17b26faa9c10f2ba563c62d2640c0ab5e4317b0d 100644
--- a/Geo/MTetrahedron.h
+++ b/Geo/MTetrahedron.h
@@ -289,8 +289,6 @@ class MTetrahedron10 : public MTetrahedron {
*
*/
-
-
/* tet order 3
* 2
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index 3d1e5a461157370e0cf3c230e9df4726bfe61f7d..09bdbecd0f95d671d4650a1ac2014ca4516f2393 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -108,12 +108,12 @@ void MVertex::writeMSH(FILE *fp, bool binary, bool saveParametric, double scalin
fprintf(fp, "\n");
}
}
+
void MVertex::writePLY2(FILE *fp)
{
if(_index < 0) return; // negative index vertices are never saved
- fprintf(fp, "%.16g %.16g %.16g\n",
- x(), y(), z());
+ fprintf(fp, "%.16g %.16g %.16g\n", x(), y(), z());
}
void MVertex::writeVRML(FILE *fp, double scalingFactor)
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index 7a6ea597a99faf3ac191465bc7ec59270f012eb9..d09a3faed6203542bcc388b92c2d11504ac01976 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -297,7 +297,6 @@ static fullMatrix<double> generatePascalPrism(int order)
return monomials;
}
-
static int nbdoftriangle(int order) { return (order + 1) * (order + 2) / 2; }
//static int nbdoftriangleserendip(int order) { return 3 * order; }
@@ -706,8 +705,6 @@ static fullMatrix<double> gmshGeneratePointsPrism(int order, bool serendip)
}
}
-// point.print("Pri ipts");
-
return point;
}
@@ -1333,19 +1330,18 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsQuad(10,true);
generate2dEdgeClosure(F.closures, 10, 4);
break;
- case MSH_PRI_6 : // first order
+ case MSH_PRI_6 :
F.numFaces = 5;
F.monomials = generatePascalPrism(1);
F.points = gmshGeneratePointsPrism(1, false);
generate3dFaceClosurePrism(F.closures, 1);
break;
- case MSH_PRI_18 : // second order
+ case MSH_PRI_18 :
F.numFaces = 5;
F.monomials = generatePascalPrism(2);
F.points = gmshGeneratePointsPrism(2, false);
generate3dFaceClosurePrism(F.closures, 2);
break;
-
default :
Msg::Error("Unknown function space %d: reverting to TET_4", tag);
F.numFaces = 4;
@@ -1357,13 +1353,14 @@ const polynomialBasis *polynomialBases::find(int tag)
F.type = tag;
F.coefficients = generateLagrangeMonomialCoefficients(F.monomials, F.points);
-// printf("Case: %d coeffs:\n",tag);
-// for (int i = 0; i<F.coefficients.size1(); i++) {
-// for (int j = 0; j<F.coefficients.size2(); j++) {
-// printf("%4.1f ",F.coefficients(i,j));
-// }
-// printf("\n");
-// }
+
+ // printf("Case: %d coeffs:\n",tag);
+ // for (int i = 0; i<F.coefficients.size1(); i++) {
+ // for (int j = 0; j<F.coefficients.size2(); j++) {
+ // printf("%4.1f ",F.coefficients(i,j));
+ // }
+ // printf("\n");
+ // }
fs.insert(std::make_pair(tag, F));
return &fs[tag];
@@ -1396,7 +1393,8 @@ const fullMatrix<double> &polynomialBases::findInjector(int tag1, int tag2)
}
#include "Bindings.h"
-void polynomialBasis::registerBindings(binding *b) {
+void polynomialBasis::registerBindings(binding *b)
+{
classBinding *cb = b->addClass<polynomialBasis>("polynomialBasis");
cb->setDescription("polynomial shape functions for elements");
methodBinding *mb = cb->addMethod
diff --git a/Solver/dofManager.h b/Solver/dofManager.h
index cb0b5fa8d9c8374f346917a6aeccfbcd77a5924d..baed02f0e0963d2723b4fd117fe17c46646514f6 100644
--- a/Solver/dofManager.h
+++ b/Solver/dofManager.h
@@ -11,14 +11,14 @@
#include <complex>
#include <map>
#include <list>
+#include <iostream>
#include "GmshConfig.h"
#include "MVertex.h"
#include "linearSystem.h"
#include "fullMatrix.h"
-#include <iostream>
-#ifdef HAVE_MPI
- #include "mpi.h"
+#if defined(HAVE_MPI)
+#include "mpi.h"
#endif
class Dof{
@@ -39,7 +39,6 @@ class Dof{
i1 = t % 10000;
i2 = t / 10000;
}
-
bool operator < (const Dof &other) const
{
if(_entity < other._entity) return true;
@@ -56,15 +55,24 @@ template<class T> struct dofTraits
{
typedef T VecType;
typedef T MatType;
- inline static void gemm (VecType &r, const MatType &m, const VecType &v, double alpha, double beta) { r = beta*r+alpha*m*v;}
+ inline static void gemm(VecType &r, const MatType &m, const VecType &v,
+ double alpha, double beta)
+ {
+ r = beta * r + alpha * m * v;
+ }
};
template<class T> struct dofTraits<fullMatrix<T> >
{
typedef fullMatrix<T> VecType;
typedef fullMatrix<T> MatType;
- inline static void gemm (VecType &r, const MatType &m, const VecType &v, double alpha, double beta) { r.gemm(m, v, alpha,beta);}
+ inline static void gemm (VecType &r, const MatType &m, const VecType &v,
+ double alpha, double beta)
+ {
+ r.gemm(m, v, alpha,beta);
+ }
};
+
/*
template<> struct dofTraits<fullVector<std::complex<double> > >
{
@@ -72,6 +80,7 @@ template<> struct dofTraits<fullVector<std::complex<double> > >
typedef fullMatrix<std::complex<double> > MatType;
};
*/
+
template<class T>
class DofAffineConstraint{
public:
@@ -111,8 +120,8 @@ class dofManager{
std::map<const std::string, linearSystem<dataMat>*> _linearSystems;
linearSystem<dataMat> *_current;
- // ********** parallel section
- private :
+ // parallel section
+ private :
// those dof are images of ghost located on another proc (id givent by the map).
// this is a first try, maybe not the final implementation
std::map<Dof, int > ghosted; // dof => procId
@@ -121,14 +130,18 @@ class dofManager{
bool _isParallel;
public:
void _parallelFinalize();
- // **********
+
public:
- dofManager(linearSystem<dataMat> *l, bool isParallel=false) : _current(l),_isParallel(isParallel), _parallelFinalized(false) { _linearSystems["A"] = l; }
- dofManager(linearSystem<dataMat> *l1, linearSystem<dataMat> *l2) : _current(l1),_isParallel(false), _parallelFinalized(false) {
+ dofManager(linearSystem<dataMat> *l, bool isParallel=false)
+ : _current(l), _isParallel(isParallel), _parallelFinalized(false)
+ {
+ _linearSystems["A"] = l;
+ }
+ dofManager(linearSystem<dataMat> *l1, linearSystem<dataMat> *l2)
+ : _current(l1), _isParallel(false), _parallelFinalized(false)
+ {
_linearSystems.insert(std::make_pair("A", l1));
_linearSystems.insert(std::make_pair("B", l2));
- //_linearSystems.["A"] = l1;
- //_linearSystems["B"] = l2;
}
virtual inline void fixDof(Dof key, const dataVec &value)
{
@@ -144,10 +157,9 @@ class dofManager{
{
fixDof(v->getNum(), Dof::createTypeWithTwoInts(iComp, iField), value);
}
-
- inline bool isFixed(Dof key) const {
- if(fixed.find(key) != fixed.end())
- {
+ inline bool isFixed(Dof key) const
+ {
+ if(fixed.find(key) != fixed.end()){
return true;
}
return false;
@@ -160,14 +172,12 @@ class dofManager{
{
return isFixed(v->getNum(), Dof::createTypeWithTwoInts(iComp, iField));
}
-
inline void numberGhostDof (Dof key, int procId) {
if (fixed.find(key) != fixed.end()) return;
if (constraints.find(key) != constraints.end()) return;
if (ghosted.find(key) != ghosted.end()) return;
ghosted[key] = procId;
}
-
inline void numberDof(Dof key)
{
if (fixed.find(key) != fixed.end()) return;
@@ -188,7 +198,6 @@ class dofManager{
{
numberDof(v->getNum(), Dof::createTypeWithTwoInts(iComp, iField));
}
-
virtual inline void getDofValue(std::vector<Dof> &keys,std::vector<dataVec> &Vals)
{
int ndofs=keys.size();
@@ -196,7 +205,6 @@ class dofManager{
Vals.resize(originalSize+ndofs);
for (int i=0;i<ndofs;++i) getDofValue(keys[i], Vals[originalSize+i]);
}
-
virtual inline void getDofValue(Dof key, dataVec &val) const
{
{
@@ -214,22 +222,21 @@ class dofManager{
}
}
{
- typename std::map<Dof, DofAffineConstraint< dataVec > >::const_iterator it = constraints.find(key);
- if (it != constraints.end())
- {
- dataVec tmp(val);
- val = it->second.shift;
- for (unsigned i=0;i<(it->second).linear.size();i++)
- {
- std::map<Dof, int>::const_iterator itu = unknown.find(((it->second).linear[i]).first);
- getDofValue(((it->second).linear[i]).first, tmp);
- dofTraits<T>::gemm(val,((it->second).linear[i]).second, tmp, 1, 1);
- }
- return ;
+ typename std::map<Dof, DofAffineConstraint< dataVec > >::const_iterator it =
+ constraints.find(key);
+ if (it != constraints.end()){
+ dataVec tmp(val);
+ val = it->second.shift;
+ for (unsigned i=0;i<(it->second).linear.size();i++){
+ std::map<Dof, int>::const_iterator itu = unknown.find
+ (((it->second).linear[i]).first);
+ getDofValue(((it->second).linear[i]).first, tmp);
+ dofTraits<T>::gemm(val,((it->second).linear[i]).second, tmp, 1, 1);
+ }
+ return ;
}
}
}
-
inline void getDofValue(int ent, int type, dataVec &v) const
{
Dof key(ent, type);
@@ -247,28 +254,26 @@ class dofManager{
return;
}
}
- {
- typename std::map<Dof, DofAffineConstraint< dataVec > >::const_iterator it = constraints.find(key);
- if (it != constraints.end())
- {
- v = it->second.shift;
- dataVec tmp(v);
- for (unsigned i=0;i<(it->second).linear.size();i++)
- {
- std::map<Dof, int>::const_iterator itu = unknown.find(((it->second).linear[i]).first);
- getDofValue(((it->second).linear[i]).first, tmp);
- dofTraits<T>::gemm(v, ((it->second).linear[i]).second, tmp, 1, 1);
- }
- return ;
+ {
+ typename std::map<Dof, DofAffineConstraint< dataVec > >::const_iterator it =
+ constraints.find(key);
+ if (it != constraints.end()){
+ v = it->second.shift;
+ dataVec tmp(v);
+ for (unsigned i=0;i<(it->second).linear.size();i++){
+ std::map<Dof, int>::const_iterator itu = unknown.find
+ (((it->second).linear[i]).first);
+ getDofValue(((it->second).linear[i]).first, tmp);
+ dofTraits<T>::gemm(v, ((it->second).linear[i]).second, tmp, 1, 1);
+ }
+ return ;
}
}
-
}
inline void getDofValue(MVertex *v, int iComp, int iField, dataVec &value) const
{
getDofValue(v->getNum(), Dof::createTypeWithTwoInts(iComp, iField), value);
}
-
inline void assemble(const Dof &R, const Dof &C, const dataMat &value)
{
if (_isParallel && !_parallelFinalized) _parallelFinalize();
@@ -278,7 +283,7 @@ class dofManager{
std::map<Dof, int>::iterator itC = unknown.find(C);
if (itC != unknown.end()){
_current->addToMatrix(itR->second, itC->second, value);
- }
+ }
else{
typename std::map<Dof, dataVec>::iterator itFixed = fixed.find(C);
if (itFixed != fixed.end()) {
@@ -294,7 +299,8 @@ class dofManager{
assembleLinConst(R, C, value);
}
}
- inline void assemble(std::vector<Dof> &R, std::vector<Dof> &C, const fullMatrix<dataMat> &m)
+ inline void assemble(std::vector<Dof> &R, std::vector<Dof> &C,
+ const fullMatrix<dataMat> &m)
{
if (_isParallel && !_parallelFinalized) _parallelFinalize();
if (!_current->isAllocated()) _current->allocate(sizeOfR());
@@ -328,74 +334,58 @@ class dofManager{
}
}
}
- else
- {
- for (unsigned int j = 0; j < C.size(); j++)
- {
- assembleLinConst(R[i], C[j], m(i, j));
+ else{
+ for (unsigned int j = 0; j < C.size(); j++){
+ assembleLinConst(R[i], C[j], m(i, j));
}
}
}
}
-
- virtual inline void assemble(std::vector<Dof> &R, const fullVector<dataMat> &m) // for linear forms
+ // for linear forms
+ virtual inline void assemble(std::vector<Dof> &R, const fullVector<dataMat> &m)
{
if (_isParallel && !_parallelFinalized) _parallelFinalize();
if (!_current->isAllocated()) _current->allocate(sizeOfR());
std::vector<int> NR(R.size());
- for (unsigned int i = 0; i < R.size(); i++)
- {
+ for (unsigned int i = 0; i < R.size(); i++){
std::map<Dof, int>::iterator itR = unknown.find(R[i]);
if (itR != unknown.end()) NR[i] = itR->second;
else NR[i] = -1;
}
- for (unsigned int i = 0; i < R.size(); i++)
- {
- if (NR[i] != -1)
- {
+ for (unsigned int i = 0; i < R.size(); i++){
+ if (NR[i] != -1){
_current->addToRightHandSide(NR[i], m(i));
}
- else
- {
- typename std::map<Dof,DofAffineConstraint<dataVec> >::iterator itConstraint;
- itConstraint = constraints.find(R[i]);
- if (itConstraint != constraints.end())
- {
- for (unsigned j=0;j<(itConstraint->second).linear.size();j++)
- {
- dataMat tmp;
- dofTraits<T>::gemm(tmp,(itConstraint->second).linear[j].second,m(i), 1, 0);
- assemble((itConstraint->second).linear[j].first,tmp);
- }
- }
- }
+ else{
+ typename std::map<Dof,DofAffineConstraint<dataVec> >::iterator itConstraint;
+ itConstraint = constraints.find(R[i]);
+ if (itConstraint != constraints.end()){
+ for (unsigned j = 0; j < (itConstraint->second).linear.size(); j++){
+ dataMat tmp;
+ dofTraits<T>::gemm(tmp,(itConstraint->second).linear[j].second,m(i), 1, 0);
+ assemble((itConstraint->second).linear[j].first,tmp);
+ }
+ }
+ }
}
}
-
-
virtual inline void assemble(std::vector<Dof> &R, const fullMatrix<dataMat> &m)
{
if (_isParallel && !_parallelFinalized) _parallelFinalize();
if (!_current->isAllocated()) _current->allocate(sizeOfR());
std::vector<int> NR(R.size());
- for (unsigned int i = 0; i < R.size(); i++)
- {
+ for (unsigned int i = 0; i < R.size(); i++){
std::map<Dof, int>::iterator itR = unknown.find(R[i]);
if (itR != unknown.end()) NR[i] = itR->second;
else NR[i] = -1;
}
- for (unsigned int i = 0; i < R.size(); i++)
- {
- if (NR[i] != -1)
- {
- for (unsigned int j = 0; j < R.size(); j++)
- {
- if (NR[j] != -1)
- {
+ for (unsigned int i = 0; i < R.size(); i++){
+ if (NR[i] != -1){
+ for (unsigned int j = 0; j < R.size(); j++){
+ if (NR[j] != -1){
_current->addToMatrix(NR[i], NR[j], m(i, j));
}
- else
- {
+ else{
typename std::map<Dof, dataVec>::iterator itFixed = fixed.find(R[j]);
if (itFixed != fixed.end()){
// tmp = -m(i,j) * itFixed->second
@@ -406,15 +396,13 @@ class dofManager{
}
}
}
- else
- {
+ else{
for (unsigned int j = 0; j < R.size(); j++){
- assembleLinConst(R[i],R[j],m(i,j));
+ assembleLinConst(R[i],R[j],m(i,j));
}
}
}
}
-
inline void assemble(int entR, int typeR, int entC, int typeC, const dataMat &value)
{
assemble(Dof(entR, typeR), Dof(entC, typeC), value);
@@ -435,20 +423,17 @@ class dofManager{
if(itR != unknown.end()){
_current->addToRightHandSide(itR->second, value);
}
- else
- {
- typename std::map<Dof,DofAffineConstraint<dataVec> >::iterator itConstraint;
- itConstraint = constraints.find(R);
- if (itConstraint != constraints.end())
- {
- for (unsigned j=0;j<(itConstraint->second).linear.size();j++)
- {
- dataMat tmp;
- dofTraits<T>::gemm(tmp,(itConstraint->second).linear[j].second,value, 1, 0);
- assemble((itConstraint->second).linear[j].first,tmp);
- }
- }
- }
+ else{
+ typename std::map<Dof,DofAffineConstraint<dataVec> >::iterator itConstraint;
+ itConstraint = constraints.find(R);
+ if (itConstraint != constraints.end()){
+ for (unsigned j = 0; j < (itConstraint->second).linear.size(); j++){
+ dataMat tmp;
+ dofTraits<T>::gemm(tmp,(itConstraint->second).linear[j].second,value, 1, 0);
+ assemble((itConstraint->second).linear[j].first,tmp);
+ }
+ }
+ }
}
inline void assemble(int entR, int typeR, const dataMat &value)
{
@@ -467,14 +452,16 @@ class dofManager{
_current->zeroMatrix();
_current->zeroRightHandSide();
}
- inline void setCurrentMatrix(std::string name){
- typename std::map<const std::string, linearSystem<dataMat>*>::iterator it = _linearSystems.find(name);
- if(it != _linearSystems.end())
- _current = it->second;
- else{
- Msg::Error("Current matrix %s not found ", name.c_str());
- throw;
- }
+ inline void setCurrentMatrix(std::string name)
+ {
+ typename std::map<const std::string, linearSystem<dataMat>*>::iterator it =
+ _linearSystems.find(name);
+ if(it != _linearSystems.end())
+ _current = it->second;
+ else{
+ Msg::Error("Current matrix %s not found ", name.c_str());
+ throw;
+ }
}
linearSystem<dataMat> *getLinearSystem(std::string &name)
{
@@ -485,13 +472,11 @@ class dofManager{
else
return 0;
}
-
inline void setLinearConstraint (Dof key, DofAffineConstraint<dataVec> &affineconstraint)
{
- constraints[key] = affineconstraint;
- // constraints.insert(std::make_pair(key,affineconstraint));
+ constraints[key] = affineconstraint;
+ // constraints.insert(std::make_pair(key,affineconstraint));
}
-
inline void assembleLinConst(const Dof &R, const Dof &C, const dataMat &value)
{
std::map<Dof, int>::iterator itR = unknown.find(R);
@@ -499,12 +484,10 @@ class dofManager{
{
typename std::map<Dof,DofAffineConstraint<dataVec> >::iterator itConstraint;
itConstraint = constraints.find(C);
- if (itConstraint != constraints.end())
- {
+ if (itConstraint != constraints.end()){
dataMat tmp(value);
- for (unsigned i=0;i<(itConstraint->second).linear.size();i++)
- {
- dofTraits<T>::gemm(tmp,(itConstraint->second).linear[i].second,value, 1, 0);
+ for (unsigned i = 0; i < (itConstraint->second).linear.size(); i++){
+ dofTraits<T>::gemm(tmp, (itConstraint->second).linear[i].second, value, 1, 0);
assemble(R,(itConstraint->second).linear[i].first,tmp);
}
dataMat tmp2(value);
@@ -512,33 +495,32 @@ class dofManager{
_current->addToRightHandSide(itR->second, tmp2);
}
}
- else // test function ; (no shift ?)
- {
+ else{ // test function ; (no shift ?)
typename std::map<Dof,DofAffineConstraint<dataVec> >::iterator itConstraint;
itConstraint = constraints.find(R);
- if (itConstraint != constraints.end())
- {
+ if (itConstraint != constraints.end()){
dataMat tmp(value);
- for (unsigned i=0;i<(itConstraint->second).linear.size();i++)
- {
- dofTraits<T>::gemm(tmp,itConstraint->second.linear[i].second,value, 1, 0);
- assemble((itConstraint->second).linear[i].first,C,tmp);
+ for (unsigned i = 0; i < (itConstraint->second).linear.size(); i++){
+ dofTraits<T>::gemm(tmp, itConstraint->second.linear[i].second, value, 1, 0);
+ assemble((itConstraint->second).linear[i].first, C, tmp);
}
}
}
}
- void getFixedDof(std::vector<Dof> &R){
+ void getFixedDof(std::vector<Dof> &R)
+ {
R.clear();
R.reserve(fixed.size());
typename std::map<Dof, dataVec>::iterator it;
- for(it=fixed.begin(); it!=fixed.end();++it){
- R.push_back (it->first);
+ for(it = fixed.begin(); it != fixed.end(); ++it){
+ R.push_back(it->first);
}
}
};
template<class T>
-void dofManager<T>::_parallelFinalize() {
+void dofManager<T>::_parallelFinalize()
+{
#ifdef HAVE_MPI
int _numStart;
int _numTotal;
@@ -546,7 +528,8 @@ void dofManager<T>::_parallelFinalize() {
_localSize = unknown.size();
if (Msg::GetCommRank() == 0){
_numStart = 0;
- }else
+ }
+ else
MPI_Recv (&_numStart, 1, MPI_INT, Msg::GetCommRank()-1, 0, MPI_COMM_WORLD, &status);
_numTotal = _numStart + _localSize;
if (Msg::GetCommRank() != Msg::GetCommSize()-1)
@@ -596,15 +579,18 @@ void dofManager<T>::_parallelFinalize() {
}
}
int index;
- while (MPI_Waitany (2*Msg::GetCommSize(), reqRecv0, &index, &status) == 0 && index != MPI_UNDEFINED) {
+ while (MPI_Waitany (2*Msg::GetCommSize(), reqRecv0, &index, &status) == 0 &&
+ index != MPI_UNDEFINED) {
if (status.MPI_TAG == 0) {
for (int j = 0; j < nRequested[index]; j++) {
- std::map<Dof,int>::iterator it = unknown.find (Dof(recv0[index][j*2], recv0[index][j*2+1]));
+ std::map<Dof,int>::iterator it = unknown.find
+ (Dof(recv0[index][j*2], recv0[index][j*2+1]));
if (it == unknown.end ())
Msg::Error ("ghost Dof does not exist on parent process");
send1[index][j] = it->second;
}
- MPI_Isend (send1[index], nRequested[index], MPI_INT, index, 1, MPI_COMM_WORLD, &reqSend1[index]);
+ MPI_Isend(send1[index], nRequested[index], MPI_INT, index, 1,
+ MPI_COMM_WORLD, &reqSend1[index]);
}
}
for (int i = 0; i<Msg::GetCommSize(); i++)
diff --git a/Solver/function.h b/Solver/function.h
index 2bc228aa7590430c05fbcc15b38816efe9572939..48e6ef5adb6487c85eeeeb36b45b37544ce9153f 100644
--- a/Solver/function.h
+++ b/Solver/function.h
@@ -18,7 +18,6 @@ struct functionReplaceCache;
class MElement;
class binding;
-
// An abstract interface to functions
class function {
public: