diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index 45d42dcc57f03d9a220ef60f57320bfa9e3fad2a..cc852b39ae4e44d059abef34d8a5f22ba8d199e1 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -62,7 +62,7 @@ void PrintUsage(const char *name)
Msg::Direct(" -saveall Save all elements (discard physical group definitions)");
Msg::Direct(" -o file Specify mesh output file name");
Msg::Direct(" -format string Set output mesh format (msh, msh1, msh2, unv, vrml, stl, mesh,");
- Msg::Direct(" bdf, p3d, cgns, med)");
+ Msg::Direct(" bdf, p3d, cgns, med, fea)");
Msg::Direct(" -bin Use binary format when available");
Msg::Direct(" -parametric Save vertices with their parametric coordinates");
Msg::Direct(" -numsubedges Set the number of subdivisions when displaying high order elements");
@@ -451,7 +451,9 @@ void GetOptions(int argc, char *argv[])
CTX::instance()->mesh.format = FORMAT_DIFF;
else if(!strcmp(argv[i], "med"))
CTX::instance()->mesh.format = FORMAT_MED;
- else
+ else if(!strcmp(argv[i], "fea"))
+ CTX::instance()->mesh.format = FORMAT_VISUALFEA;
+ else
Msg::Fatal("Unknown mesh format");
i++;
}
diff --git a/Common/CreateFile.cpp b/Common/CreateFile.cpp
index 541375f735dbb0b752dc6f06477d9faccb04c7c5..7d369001c5dad8bed5a0548d82a0a7d906fca1d5 100644
--- a/Common/CreateFile.cpp
+++ b/Common/CreateFile.cpp
@@ -38,6 +38,7 @@ int GuessFileFormatFromFileName(std::string fileName)
else if(ext == ".stl") return FORMAT_STL;
else if(ext == ".cgns") return FORMAT_CGNS;
else if(ext == ".med") return FORMAT_MED;
+ else if(ext == ".fea") return FORMAT_VISUALFEA;
else if(ext == ".mesh") return FORMAT_MESH;
else if(ext == ".bdf") return FORMAT_BDF;
else if(ext == ".diff") return FORMAT_DIFF;
@@ -78,6 +79,7 @@ std::string GetDefaultFileName(int format)
case FORMAT_STL: name += ".stl"; break;
case FORMAT_CGNS: name += ".cgns"; break;
case FORMAT_MED: name += ".med"; break;
+ case FORMAT_VISUALFEA: name += ".fea"; break;
case FORMAT_MESH: name += ".mesh"; break;
case FORMAT_BDF: name += ".bdf"; break;
case FORMAT_DIFF: name += ".diff"; break;
@@ -145,6 +147,7 @@ static PixelBuffer *GetCompositePixelBuffer(GLenum format, GLenum type)
void CreateOutputFile(std::string fileName, int format)
{
+
if(fileName.empty())
fileName = GetDefaultFileName(format);
@@ -204,6 +207,12 @@ void CreateOutputFile(std::string fileName, int format)
CTX::instance()->mesh.saveAll, CTX::instance()->mesh.scalingFactor);
break;
+ case FORMAT_VISUALFEA:
+ GModel::current()->writeVisualFEA
+ (fileName, CTX::instance()->mesh.saveElementTagType,
+ CTX::instance()->mesh.saveAll, CTX::instance()->mesh.scalingFactor);
+ break;
+
case FORMAT_BDF:
GModel::current()->writeBDF
(fileName, CTX::instance()->mesh.bdfFieldFormat,
diff --git a/Common/GmshDefines.h b/Common/GmshDefines.h
index fa4bc565c6df4bc5e0f47433175ab4d2355ac182..d4b66d6be42388a5195565d8897cb9b50dbec696 100644
--- a/Common/GmshDefines.h
+++ b/Common/GmshDefines.h
@@ -40,6 +40,7 @@
#define FORMAT_BREP 35
#define FORMAT_STEP 36
#define FORMAT_IGES 37
+#define FORMAT_VISUALFEA 38
// Element types
#define TYPE_PNT 1
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index fda0b1799e0c0a3080f1149ee98f58fdf817d198..fc3b57ad2ce1c62806db8133afa7cd3988ea6777 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -23,6 +23,7 @@
#include "crossConfTerm.h"
#include "convexCombinationTerm.h"
#include "linearSystemGMM.h"
+#include "linearSystemTAUCS.h"
#include "linearSystemCSR.h"
#include "linearSystemFull.h"
@@ -481,12 +482,9 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
if (!_lsys) {
#if defined(HAVE_GMM)
- linearSystemGmm<double> *_lsysb = new linearSystemGmm<double>;
- //linearSystemCSRGmm<double> lsys;
- _lsysb->setPrec(1.e-15);
+ linearSystemGmm<double> * _lsysb = new linearSystemGmm<double>;
_lsysb->setGmres(1);
- if(Msg::GetVerbosity() == 99) _lsysb->setNoisy(2);
- _lsys = _lsysb;
+ _lsys=_lsysb;
#else
_lsys = new linearSystemFull<double>;
#endif
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 9aab297c182972cf498d34627837d6a5564da9a8..c7fffad427bcc54ae19cbf31319dd0f070590ea6 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -344,6 +344,10 @@ class GModel
int writeMSH(const std::string &name, double version=1.0, bool binary=false,
bool saveAll=false, bool saveParametric=false, double scalingFactor=1.0);
+ // VisualFEA file format
+ int writeVisualFEA(const std::string &name, int elementTagType,
+ bool saveAll, double scalingFactor);
+
// mesh statistics (saved as a Gmsh post-processing view)
int writePOS(const std::string &name, bool printElementary,
bool printElementNumber, bool printGamma, bool printEta, bool printRho,
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 8cb5618185b88ad8ec46be7e52148575bc2d5f2f..a9a7fd8adbab390e2d73cf0a335026f6fe9efd43 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -1606,6 +1606,97 @@ int GModel::readMESH(const std::string &name)
return 1;
}
+
+int GModel::writeVisualFEA(const std::string &name, int elementTagType,
+ bool saveAll, double scalingFactor)
+{
+ FILE *fp = fopen(name.c_str(), "w");
+ if(!fp){
+ Msg::Error("Unable to open file '%s'", name.c_str());
+ return 0;
+ }
+
+ if(noPhysicalGroups()) saveAll = true;
+
+ int numVertices = indexMeshVertices(saveAll);
+ int numHexahedra = 0, numTriangles = 0, numQuadrangles = 0, numTetrahedra = 0;
+ for(fiter it = firstFace(); it != lastFace(); ++it){
+ if(saveAll || (*it)->physicals.size()){
+ numTriangles += (*it)->triangles.size();
+ numQuadrangles += (*it)->quadrangles.size();
+ }
+ }
+ for(riter it = firstRegion(); it != lastRegion(); ++it){
+ if(saveAll || (*it)->physicals.size()){
+ numTetrahedra += (*it)->tetrahedra.size();
+ numHexahedra += (*it)->hexahedra.size();
+ }
+ }
+
+ fprintf(fp,"33\n");
+ fprintf(fp,"%d %d %d\n", numVertices, numTriangles+numQuadrangles, numTetrahedra+numHexahedra);
+
+ std::vector<GEntity*> entities;
+ getEntities(entities);
+ for(unsigned int i = 0; i < entities.size(); i++)
+ for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++)
+ fprintf(fp,"%d %g %g %g\n",entities[i]->mesh_vertices[j]->getNum(),
+ entities[i]->mesh_vertices[j]->x()*scalingFactor,
+ entities[i]->mesh_vertices[j]->y()*scalingFactor,
+ entities[i]->mesh_vertices[j]->z()*scalingFactor);
+ int iElement = 0;
+ for(fiter it = firstFace(); it != lastFace(); ++it){
+ int numPhys = (*it)->physicals.size();
+ if(saveAll || numPhys){
+ for(unsigned int i = 0; i < (*it)->triangles.size(); i++){
+ MTriangle *t = (*it)->triangles[i];
+ fprintf(fp,"%d %d 3 %d %d %d\n",iElement++,(*it)->tag(),
+ t->getVertex(0)->getNum(),
+ t->getVertex(1)->getNum(),
+ t->getVertex(2)->getNum());
+ }
+ for(unsigned int i = 0; i < (*it)->quadrangles.size(); i++){
+ MQuadrangle *t = (*it)->quadrangles[i];
+ fprintf(fp,"%d %d 4 %d %d %d\n",iElement++,(*it)->tag(),
+ t->getVertex(0)->getNum(),
+ t->getVertex(1)->getNum(),
+ t->getVertex(2)->getNum(),
+ t->getVertex(3)->getNum());
+ }
+ }
+ }
+
+ for(riter it = firstRegion(); it != lastRegion(); ++it){
+ int numPhys = (*it)->physicals.size();
+ if(saveAll || numPhys){
+ for(unsigned int i = 0; i < (*it)->tetrahedra.size(); i++){
+ MTetrahedron *t = (*it)->tetrahedra[i];
+ fprintf(fp,"%d %d 3 %d %d %d\n",iElement++,(*it)->tag(),
+ t->getVertex(0)->getNum(),
+ t->getVertex(1)->getNum(),
+ t->getVertex(2)->getNum(),
+ t->getVertex(3)->getNum());
+ }
+ for(unsigned int i = 0; i < (*it)->hexahedra.size(); i++){
+ MHexahedron *t = (*it)->hexahedra[i];
+ fprintf(fp,"%d %d 3 %d %d %d\n",iElement++,(*it)->tag(),
+ t->getVertex(0)->getNum(),
+ t->getVertex(1)->getNum(),
+ t->getVertex(2)->getNum(),
+ t->getVertex(3)->getNum(),
+ t->getVertex(4)->getNum(),
+ t->getVertex(5)->getNum(),
+ t->getVertex(6)->getNum(),
+ t->getVertex(7)->getNum());
+ }
+ }
+ }
+
+ fclose(fp);
+ return 1;
+}
+
+
int GModel::writeMESH(const std::string &name, int elementTagType,
bool saveAll, double scalingFactor)
{
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 54a76fc919b08886d8929dab2e82d7f9fb4a741f..8407e8b1f52ed866b1e76ac5c9839a79a19b4b7e 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -106,6 +106,9 @@ class MElement
// get the edges
virtual int getNumEdges() = 0;
virtual MEdge getEdge(int num) = 0;
+ // give an MEdge as input and get its local number and sign
+ virtual void getEdgeInfo (const MEdge & edge, int &ithEdge, int &sign) const
+ {throw;}
// get an edge representation for drawing
virtual int getNumEdgesRep() = 0;
@@ -120,6 +123,9 @@ class MElement
// get the faces
virtual int getNumFaces() = 0;
virtual MFace getFace(int num) = 0;
+ // give an MFace as input and get its local number, sign and rotation
+ virtual void getFaceInfo (const MFace & face, int &ithFace, int &sign, int &rot)const
+ {throw;}
// get a face representation for drawing
virtual int getNumFacesRep() = 0;
diff --git a/Geo/MQuadrangle.h b/Geo/MQuadrangle.h
index 7d009b041bdb258291a885f81d8b41479c34f0f3..376265d3fbface8cb53dbb095d258bdf651eab56 100644
--- a/Geo/MQuadrangle.h
+++ b/Geo/MQuadrangle.h
@@ -68,6 +68,19 @@ class MQuadrangle : public MElement {
{
return MEdge(_v[edges_quad(num, 0)], _v[edges_quad(num, 1)]);
}
+ virtual void getEdgeInfo (const MEdge & edge, int &ithEdge, int &sign) const {
+ for (ithEdge=0;ithEdge<4;ithEdge++){
+ const MVertex *v0 = _v[edges_quad(ithEdge, 0)];
+ const MVertex *v1 = _v[edges_quad(ithEdge, 1)];
+ if (v0 == edge.getVertex(0) && v1 == edge.getVertex(1)){
+ sign = 1; return;
+ }
+ if (v1 == edge.getVertex(0) && v0 == edge.getVertex(1)){
+ sign = -1; return;
+ }
+ }
+ throw;
+ }
virtual int getNumEdgesRep(){ return 4; }
virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
{
diff --git a/Geo/MTetrahedron.cpp b/Geo/MTetrahedron.cpp
index f8b51640869a9eedce205cd31a8224a91b268f1d..cac4b488a559ad9c54c6f57e0c771c97d4bb7e06 100644
--- a/Geo/MTetrahedron.cpp
+++ b/Geo/MTetrahedron.cpp
@@ -231,3 +231,30 @@ void MTetrahedron::getIntegrationPoints(int pOrder, int *npts, IntPt **pts) cons
*npts = getNGQTetPts(pOrder);
*pts = getGQTetPts(pOrder);
}
+void MTetrahedron::getFaceInfo (const MFace & face, int &ithFace, int &sign, int &rot)const{
+ for (ithFace=0;ithFace<4;ithFace++){
+ MVertex *v0 = _v[faces_tetra(ithFace,0)];
+ MVertex *v1 = _v[faces_tetra(ithFace,1)];
+ MVertex *v2 = _v[faces_tetra(ithFace,2)];
+
+ if (v0 == face.getVertex(0) && v1 == face.getVertex(1) && v2 == face.getVertex(2)){
+ sign = 1; rot = 0; return;
+ }
+ if (v0 == face.getVertex(1) && v1 == face.getVertex(2) && v2 == face.getVertex(0)){
+ sign = 1; rot = 1; return;
+ }
+ if (v0 == face.getVertex(2) && v1 == face.getVertex(0) && v2 == face.getVertex(1)){
+ sign = 1; rot = 2; return;
+ }
+ if (v0 == face.getVertex(0) && v1 == face.getVertex(2) && v2 == face.getVertex(1)){
+ sign = -1; rot = 0; return;
+ }
+ if (v0 == face.getVertex(1) && v1 == face.getVertex(0) && v2 == face.getVertex(2)){
+ sign = -1; rot = 1; return;
+ }
+ if (v0 == face.getVertex(2) && v1 == face.getVertex(1) && v2 == face.getVertex(0)){
+ sign = -1; rot = 2; return;
+ }
+ }
+ throw;
+}
diff --git a/Geo/MTetrahedron.h b/Geo/MTetrahedron.h
index cf15ea0224882a3b4930dcffce604036fe58ecbc..0d0c2745071069a4e935e47a44798cc29eea6085 100644
--- a/Geo/MTetrahedron.h
+++ b/Geo/MTetrahedron.h
@@ -89,6 +89,7 @@ class MTetrahedron : public MElement {
_v[faces_tetra(num, 1)],
_v[faces_tetra(num, 2)]);
}
+ virtual void getFaceInfo (const MFace & face, int &ithFace, int &sign, int &rot)const;
virtual int getNumFacesRep(){ return 4; }
virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
{
diff --git a/Geo/MTriangle.h b/Geo/MTriangle.h
index fef5d9d02e915856bbbd02e7d9db913e08e56d81..c63c328c04f89a37d797f73b86bac785951a19df 100644
--- a/Geo/MTriangle.h
+++ b/Geo/MTriangle.h
@@ -71,6 +71,19 @@ class MTriangle : public MElement {
{
return MEdge(_v[edges_tri(num, 0)], _v[edges_tri(num, 1)]);
}
+ virtual void getEdgeInfo (const MEdge & edge, int &ithEdge, int &sign) const {
+ for (ithEdge=0;ithEdge<3;ithEdge++){
+ const MVertex *v0 = _v[edges_tri(ithEdge, 0)];
+ const MVertex *v1 = _v[edges_tri(ithEdge, 1)];
+ if (v0 == edge.getVertex(0) && v1 == edge.getVertex(1)){
+ sign = 1; return;
+ }
+ if (v1 == edge.getVertex(0) && v0 == edge.getVertex(1)){
+ sign = -1; return;
+ }
+ }
+ throw;
+ }
virtual int getNumEdgesRep(){ return 3; }
virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
{
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 6a2b5dbe2b7d639f0950209362089855ef4bc5da..d1a4f15ae2d7920a4492996c27c54352e99354ab 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -870,8 +870,10 @@ void bowyerWatsonFrontal(GFace *gf)
// }
}
- // char name[245];
- // sprintf(name,"frontal%d.pos", gf->tag());
- // _printTris (name, AllTris, Us, Vs);
+ char name[245];
+ sprintf(name,"frontal%d-real.pos", gf->tag());
+ _printTris (name, AllTris, Us, Vs,false);
+ sprintf(name,"frontal%d-param.pos", gf->tag());
+ _printTris (name, AllTris, Us, Vs,true);
transferDataStructure(gf, AllTris, Us, Vs);
}
diff --git a/Numeric/functionSpace.cpp b/Numeric/functionSpace.cpp
index 60dbf13ce84f542c02218038f5792e7e846d2723..fe26349a3df59284d45a7c39f95a0947d724a35c 100644
--- a/Numeric/functionSpace.cpp
+++ b/Numeric/functionSpace.cpp
@@ -542,6 +542,102 @@ static fullMatrix<double> generateLagrangeMonomialCoefficients
return coefficient;
}
+
+// -----------------------------------------------------------------------------
+static void getFaceClosure(int iFace, int iSign, int iRotate , std::vector<int> & closure, int order) {
+ closure.clear();
+ closure.reserve((order+1)*(order+2)/2);
+ switch (order){
+ case 0:
+ closure[0]=0;
+ break;
+ default:
+ int face[4][3] = {{0, 1, 2},{0, 4, -3},{1, 5, -4},{-2, 5, -3}};
+ int order1node[4][3] = {{0,1,2},{0,1,3},{1,2,3},{0,2,3}};
+ // int facedofstmp[200];
+ //face 0 | 0 1 2
+ //face 1 | 0 4 -3
+ //face 2 | 1 5 -4
+ //face 3 | -2 5 -3
+ // edge 0 | 4 -> 4+order-2
+ // edge 1 | 4+order-1 -> 4 + 2*(order-1)-1
+ // edge 2 | 4+2*(order-1)-> 4+ 3*(order-1)-1
+ // edge k | 4+k*(order-1) -> 4+(k+1)*(order-1)-1
+ for (int i =0;i<3;++i){
+ int k=(3+(iSign*i)+iRotate)%3;
+ closure[i]=order1node[iFace][k];
+ }
+ for (int i =0;i < 3;++i){
+ int edgenumber = iSign*face[iFace][(6+i*iSign+(-1+iSign)/2+iRotate)%3];
+ for (size_t k =0; k< (order-1); k++){
+ if (edgenumber > 0 || ((edgenumber == 0) && (iSign >0) ))
+ closure [3+i*(order-1)+k] = 4 + edgenumber*(order-1)+k;
+ else
+ closure [3+i*(order-1)+k] = 4 +(1- edgenumber)*(order-1)-1-k;
+ }
+ }
+ int fi = 3+3*(order-1);
+ int ti = 4+6*(order-1);
+ int ndofff = (order-3+2)*(order-3+1)/2;
+ ti = ti +iFace * ndofff;
+ for (size_t k=0; k<order/3;k++){
+ int orderint = order - 3 - k*3;
+ if (orderint>0){
+ for (int ci =0; ci <3 ; ci++){
+ int shift = (3+iSign*ci+iRotate)%3;
+ closure[fi+ci] = ti+shift;
+ }
+ fi= fi+3; ti= ti+3;
+ for (int l=0; l< orderint-1; l++){
+ for (int ei =0; ei<3; ei++)
+ {
+ int edgenumber = (6+ei*iSign+(-1+iSign)/2+iRotate)%3;
+ if (iSign > 0) closure[fi+ei*(orderint-1)+l] = ti + edgenumber*(orderint-1)+l;
+ else closure[fi+ei*(orderint-1)+l] = ti + (1+edgenumber)*(orderint-1)-1-l;
+ }
+ }
+ fi=fi+3*(orderint-1); ti = ti+3*(orderint -1);
+ }
+ else {
+ closure[fi] = ti ;
+ ti++;
+ fi++;
+ }
+ }
+ break;
+ }
+ return;
+}
+
+static void generate3dFaceClosure ( functionSpace::clCont & closure , int order) {
+ for (int iRotate=0;iRotate<3;iRotate++){
+ for (int iSign=1;iSign>=-1;iSign-=2){
+ for (int iFace=0;iFace<4;iFace++){
+ std::vector<int> closure_face;
+ getFaceClosure(iFace,iSign,iRotate,closure_face,order);
+ closure.push_back(closure_face);
+ }
+ }
+ }
+}
+
+static void generate2dEdgeClosure ( functionSpace::clCont & closure , int order) {
+ closure.clear();
+ // first give edge nodes of the three edges in direct order
+ int index = 3;
+ std::vector<int> c0,c1,c2;
+ for (int i = 0; i < order - 1; i++, index++) c0.push_back(index);
+ for (int i = 0; i < order - 1; i++, index++) c1.push_back(index);
+ for (int i = 0; i < order - 1; i++, index++) c2.push_back(index);
+ closure.push_back(c0);closure.push_back(c1);closure.push_back(c2);
+ // then give edge nodes in reverse order
+ std::vector<int> c3,c4,c5;
+ for (int i=c0.size()-1;i>=0;i--)c3.push_back(c0[i]);
+ for (int i=c1.size()-1;i>=0;i--)c4.push_back(c1[i]);
+ for (int i=c2.size()-1;i>=0;i--)c5.push_back(c2[i]);
+ closure.push_back(c3);closure.push_back(c4);closure.push_back(c5);
+}
+
std::map<int, functionSpace> functionSpaces::fs;
const functionSpace &functionSpaces::find(int tag)
@@ -579,67 +675,83 @@ const functionSpace &functionSpaces::find(int tag)
case MSH_TRI_3 :
F.monomials = generatePascalTriangle(1);
F.points = gmshGeneratePointsTriangle(1, false);
+ generate2dEdgeClosure (F.edgeClosure,1);
break;
case MSH_TRI_6 :
F.monomials = generatePascalTriangle(2);
F.points = gmshGeneratePointsTriangle(2, false);
+ generate2dEdgeClosure (F.edgeClosure,2);
break;
case MSH_TRI_9 :
F.monomials = generatePascalSerendipityTriangle(3);
F.points = gmshGeneratePointsTriangle(3, true);
+ generate2dEdgeClosure (F.edgeClosure,3);
break;
case MSH_TRI_10 :
F.monomials = generatePascalTriangle(3);
F.points = gmshGeneratePointsTriangle(3, false);
+ generate2dEdgeClosure (F.edgeClosure,3);
break;
case MSH_TRI_12 :
F.monomials = generatePascalSerendipityTriangle(4);
F.points = gmshGeneratePointsTriangle(4, true);
+ generate2dEdgeClosure (F.edgeClosure,4);
break;
case MSH_TRI_15 :
F.monomials = generatePascalTriangle(4);
F.points = gmshGeneratePointsTriangle(4, false);
+ generate2dEdgeClosure (F.edgeClosure,4);
break;
case MSH_TRI_15I :
F.monomials = generatePascalSerendipityTriangle(5);
F.points = gmshGeneratePointsTriangle(5, true);
+ generate2dEdgeClosure (F.edgeClosure,5);
break;
case MSH_TRI_21 :
F.monomials = generatePascalTriangle(5);
F.points = gmshGeneratePointsTriangle(5, false);
+ generate2dEdgeClosure (F.edgeClosure,5);
break;
case MSH_TET_4 :
F.monomials = generatePascalTetrahedron(1);
F.points = gmshGeneratePointsTetrahedron(1, false);
+ generate3dFaceClosure (F.faceClosure,1);
break;
case MSH_TET_10 :
F.monomials = generatePascalTetrahedron(2);
F.points = gmshGeneratePointsTetrahedron(2, false);
+ generate3dFaceClosure (F.faceClosure,2);
break;
case MSH_TET_20 :
F.monomials = generatePascalTetrahedron(3);
F.points = gmshGeneratePointsTetrahedron(3, false);
+ generate3dFaceClosure (F.faceClosure,3);
break;
case MSH_TET_35 :
F.monomials = generatePascalTetrahedron(4);
F.points = gmshGeneratePointsTetrahedron(4, false);
+ generate3dFaceClosure (F.faceClosure,4);
break;
case MSH_TET_34 :
F.monomials = generatePascalSerendipityTetrahedron(4);
F.points = gmshGeneratePointsTetrahedron(4, true);
+ generate3dFaceClosure (F.faceClosure,4);
break;
case MSH_TET_52 :
F.monomials = generatePascalSerendipityTetrahedron(5);
F.points = gmshGeneratePointsTetrahedron(5, true);
+ generate3dFaceClosure (F.faceClosure,5);
break;
case MSH_TET_56 :
F.monomials = generatePascalTetrahedron(5);
F.points = gmshGeneratePointsTetrahedron(5, false);
+ generate3dFaceClosure (F.faceClosure,5);
break;
default :
Msg::Error("Unknown function space %d: reverting to TET_4", tag);
F.monomials = generatePascalTetrahedron(1);
F.points = gmshGeneratePointsTetrahedron(1, false);
+ generate3dFaceClosure (F.faceClosure,1);
break;
}
F.coefficients = generateLagrangeMonomialCoefficients(F.monomials, F.points);
diff --git a/Numeric/functionSpace.h b/Numeric/functionSpace.h
index 9ba14b7d84ca57c12faf1a762819c1b7c3863208..8788652ff4b544ad4ba9bbd87b60b62cc6bfe4a6 100644
--- a/Numeric/functionSpace.h
+++ b/Numeric/functionSpace.h
@@ -8,13 +8,27 @@
#include <math.h>
#include <map>
+#include <vector>
#include "fullMatrix.h"
+// presently thos function spaces are only for simplices
+// should be extended to other elements like quads and hexes
struct functionSpace
{
+ typedef std::vector<std::vector<int> > clCont;
+ clCont faceClosure;
+ clCont edgeClosure;
fullMatrix<double> points;
fullMatrix<double> monomials;
fullMatrix<double> coefficients;
+ // for a given face/edge, with both a sign and a rotation,
+ // give an ordered list of nodes on this face/edge
+ std::vector<int> & getFaceClosure (int iFace, int iSign, int iRot){
+ return faceClosure[iFace+4*(iSign==1?0:1)+8*iRot];
+ }
+ inline std::vector<int> & getEdgeClosure (int iEdge, int iSign){
+ return edgeClosure[iSign == 1 ? iEdge : 3+iEdge];
+ }
inline void evaluateMonomials(double u, double v, double w, double p[]) const
{
for (int j = 0; j < monomials.size1(); j++) {
diff --git a/Solver/linearSystemCSR.cpp b/Solver/linearSystemCSR.cpp
index 52b9d8faf942ed26a937a617f389ac1551840ed1..0483916485941fed71cb461e532a87ed8028800d 100644
--- a/Solver/linearSystemCSR.cpp
+++ b/Solver/linearSystemCSR.cpp
@@ -245,12 +245,12 @@ static void _sort2_xkws(unsigned long n, double arr[], INDEX_TYPE ai[], INDEX_TY
}
template <class scalar>
-static void sortColumns(int NbLines,
- int nnz,
- INDEX_TYPE *ptr,
- INDEX_TYPE *jptr,
- INDEX_TYPE *ai,
- scalar *a)
+void sortColumns_(int NbLines,
+ int nnz,
+ INDEX_TYPE *ptr,
+ INDEX_TYPE *jptr,
+ INDEX_TYPE *ai,
+ scalar *a)
{
// replace pointers by lines
int *count = new int [NbLines];
@@ -290,7 +290,7 @@ template<>
int linearSystemCSRGmm<double>::systemSolve()
{
if (!sorted)
- sortColumns(_b->size(),
+ sortColumns_(_b->size(),
CSRList_Nbr(a_),
(INDEX_TYPE *) ptr_->array,
(INDEX_TYPE *) jptr_->array,
diff --git a/Solver/linearSystemTAUCS.cpp b/Solver/linearSystemTAUCS.cpp
index f2806251e2a9d7492a183e32aedb261c2b9eb1be..e84ab899fcbd800a31138d51df90fd06c38c4fa8 100644
--- a/Solver/linearSystemTAUCS.cpp
+++ b/Solver/linearSystemTAUCS.cpp
@@ -14,7 +14,7 @@ extern "C" {
}
template <class scalar>
-void sortColumns(int NbLines,
+void sortColumns_(int NbLines,
int nnz,
INDEX_TYPE *ptr,
INDEX_TYPE *jptr,
@@ -25,7 +25,7 @@ template<>
int linearSystemCSRTaucs<double>::systemSolve()
{
if(!sorted){
- sortColumns(_b->size(),
+ sortColumns_(_b->size(),
CSRList_Nbr(a_),
(INDEX_TYPE *) ptr_->array,
(INDEX_TYPE *) jptr_->array,
diff --git a/contrib/DiscreteIntegration/Integration3D.h b/contrib/DiscreteIntegration/Integration3D.h
index e993da8e8f3f3d615280b1aea83097ea98df8d8f..8cac2681b26bf35a95b2feeaf700a4d2c8084c08 100644
--- a/contrib/DiscreteIntegration/Integration3D.h
+++ b/contrib/DiscreteIntegration/Integration3D.h
@@ -90,7 +90,7 @@ class DI_Point
// --------------------------------------------------------------------------------------------------
class DI_IntegrationPoint
{
- // coordiantes
+ // coordinates
double x_, y_, z_;
// local coordinates
double xl_, yl_, zl_;