diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index 7c616f22b69812672f1dfa1eea18a5d270957636..e998901bc093fbeefd59ba0c920324e2972fefe7 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -1204,7 +1204,7 @@ StringXNumber MeshOptions_Number[] = {
{ F|O, "Recombine3DLevel" , opt_mesh_recombine3d_level , 0 ,
"3d recombination level (0: hex, 1: hex+prisms, 2: hex+prism+pyramids)" },
{ F|O, "Recombine3DConformity" , opt_mesh_recombine3d_conformity , 0 ,
- "3d recombination conformity type (0: nonconforming, 1: conformity using trihedra, 2: conformity using pyramids)" },
+ "3d recombination conformity type (0: nonconforming, 1: trihedra, 2: pyramids+trihedra, 3:pyramids+hexSplit+trihedra, 4:hexSplit+trihedra)" },
{ F|O, "DoRecombinationTest" , opt_mesh_do_recombination_test , 0 ,
"Apply recombination algorithm for test" },
{ F|O, "RecombinationTestHorizStart" , opt_mesh_recombination_test_start , 1 ,
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 743d664f774702d271e3026644855a43be8e7606..e5a1e4c254ee4edcdeef76caf51b4d2692e7243a 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -600,53 +600,16 @@ static void addToMap(std::multimap< MFace , MElement *, Less_Face> &faceToElemen
MFace outFace = fit->first;
std::vector<std::pair<MElement *, bool> > &neigh = elToNeighbors[el];
for (size_t iN = 0; iN < neigh.size(); iN++)
- if (neigh[iN].first == fit->second)
+ if (neigh[iN].first == fit->second) //Neigbor is already in the map
return;
int i0 = -1;
while (face.getVertex(0) != outFace.getVertex(++i0));
bool sameOrientation = face.getVertex(1) == outFace.getVertex((i0+1)%face.getNumVertices());
- // if (sameOrientation) printf("Non-matching orientation between el %i and el %i for face IN %i %i %i out %i %i %i\n",el->getNum(),fit->second->getNum(),face.getVertex(0)->getNum(),face.getVertex(1)->getNum(),face.getVertex(2)->getNum(),outFace.getVertex(0)->getNum(),outFace.getVertex(1)->getNum(),outFace.getVertex(2)->getNum());
neigh.push_back(std::make_pair(fit->second, !sameOrientation));
elToNeighbors[fit->second].push_back(std::make_pair(el, !sameOrientation));
}
}
-//static void checkConformity(std::multimap< MFace , MElement *, Less_Face> &faceToElement, std::map< MElement *, std::vector < std::pair <MElement *, bool> > > &elToNeighbors, const MFace &face, MElement *el){
-// //Check conformity
-// if (face.getNumVertices() == 4){
-// MVertex *v0 = face.getVertex(0);
-// MVertex *v1 = face.getVertex(1);
-// MVertex *v2 = face.getVertex(2);
-// MVertex *v3 = face.getVertex(3);
-// MFace f[5] = {MFace(v0, v1, v2), MFace(v0, v2, v3), MFace(v0, v1, v3), MFace(v1, v2, v3)};
-// if (faceToElement.find(f[0])!= faceToElement.end() || faceToElement.find(f[1])!= faceToElement.end() || faceToElement.find(f[2])!= faceToElement.end() || faceToElement.find(f[3])!= faceToElement.end() ){
-// //Check if the element is a trihedron used to recover conformity:
-// //A single element contains the two triangles, the quad and no other face
-// std::multiset<size_t> facesNeighborsId;
-// int nFaces=0;
-// for (int iFace = 0; iFace<4; iFace++){
-// int nbCount = faceToElement.count(f[iFace]);
-// if (nbCount != 0){
-// nFaces++;
-// if (nFaces > 3) Msg::Fatal("Topological fault: A quad face shares more than 3 triangles (element %i)",el->getNum());
-// std::pair <std::multimap< MFace , MElement *, Less_Face>::iterator, std::multimap< MFace , MElement *, Less_Face>::iterator> faceNeighbors;
-// faceNeighbors = faceToElement.equal_range(f[iFace]);
-// for (std::multimap< MFace , MElement *, Less_Face>::iterator itEl=faceNeighbors.first; itEl!=faceNeighbors.second; ++itEl)
-// facesNeighborsId.insert(itEl->second->getNum());
-// }
-// }
-// //In some cases, there can be 3 triangles matching a quad face and result in a valid trihedron (1 or 4 triangles not allowed)
-// if ((nFaces != 2 || nFaces != 3) && faceToElement.count(face)!= 2) Msg::Fatal("Nonconforming element which is not a trihedron (element %i nFaces %i nodes %i %i %i %i)",el->getNum(),nFaces,v0->getNum(),v1->getNum(),v2->getNum(),v3->getNum());
-// size_t maxOccurence=0;
-// for (std::multiset<size_t>::iterator iElId = facesNeighborsId.begin(); iElId != facesNeighborsId.end(); ++iElId){
-// maxOccurence = std::max(facesNeighborsId.count(*iElId), maxOccurence);
-// }
-// if (maxOccurence != 3) //A trihedron should have 3 and only 3 neighbors
-// Msg::Fatal("Nonconforming face for element: %i vertices %i %i %i %i. Number of neighbors: %i",el->getNum(), v0->getNum(), v1->getNum(), v2->getNum(), v3->getNum(),maxOccurence);
-// }
-// }
-//}
-
static void checkConformity(std::multimap< MFace , MElement *, Less_Face> &faceToElement, std::map< MElement *, std::vector < std::pair <MElement *, bool> > > &elToNeighbors, MFace face, MElement *el){
int connectivity = faceToElement.count(face);
if (ElementType::ParentTypeFromTag(el->getType()) == TYPE_TRIH){
@@ -658,8 +621,8 @@ static void checkConformity(std::multimap< MFace , MElement *, Less_Face> &faceT
for (int iV = 0; iV < face.getNumVertices(); iV++)
if (face.getVertex(iV)->onWhat()->dim() == 3 || connectivity != 1){
for (int jV = 0; jV < face.getNumVertices(); jV++)
- Msg::Info("Vertex %i",face.getVertex(jV)->getNum());
- Msg::Fatal("Non conforming element %i (nb connections for a face %i vertex %i)",el->getNum(), connectivity, face.getVertex(iV)->getNum());
+ Msg::Info("Vertex %i dim %i",face.getVertex(jV)->getNum(), face.getVertex(iV)->onWhat()->dim());
+ Msg::Fatal("Non conforming element %i (%i connection(s) for a face located on dim %i (vertex %i))",el->getNum(), connectivity,face.getVertex(iV)->onWhat()->dim(), face.getVertex(iV)->getNum());
}
}
}
diff --git a/Mesh/yamakawa.cpp b/Mesh/yamakawa.cpp
index f78f3a6bff9f6a0bf53697683613e2415b497e57..8fd30452cf721438c8ab7645e972bb69f8d77949 100644
--- a/Mesh/yamakawa.cpp
+++ b/Mesh/yamakawa.cpp
@@ -4926,11 +4926,6 @@ void PostOp::execute(GRegion* gr,int level, int conformity){
estimate2 = 0;
iterations = 0;
-
- nbFAILVERTEX=0;
- nbFAILVALIDITY=0;
-
-
build_tuples(gr);
if (level >= 2){
@@ -4940,25 +4935,35 @@ void PostOp::execute(GRegion* gr,int level, int conformity){
rearrange(gr);
}
- printf("nb FAILS VERTEX after pyramids1 %i\n",nbFAILVERTEX);
- printf("nb FAILS VALIDITY after pyramids1 %i\n",nbFAILVALIDITY);
- nbFAILVERTEX=0;
- nbFAILVALIDITY=0;
+if(conformity == 2 || conformity == 3){
+ init_markings(gr);
+ build_vertex_to_tetrahedra(gr);
+ build_vertex_to_pyramids(gr);
+ pyramids2(gr);
+ rearrange(gr);
+}
+
+ if (conformity == 3 || conformity == 4){
+ init_markings_hex(gr);
+ build_vertex_to_tetrahedra(gr);
+ build_vertex_to_pyramids(gr);
+ split_hexahedra(gr);
+ rearrange(gr);
- if(conformity == 2 || level >= 3){
- init_markings(gr);
+ init_markings_pri(gr);
build_vertex_to_tetrahedra(gr);
build_vertex_to_pyramids(gr);
- pyramids2(gr);
+ split_prisms(gr);
rearrange(gr);
- }
-
-
- printf("nb FAILS VERTEX after pyramids2 %i\n",nbFAILVERTEX);
- printf("nb FAILS VALIDITY after pyramids2 %i\n",nbFAILVALIDITY);
-
+
+ init_markings_pyr(gr);
+ build_vertex_to_tetrahedra(gr);
+ build_vertex_to_pyramids(gr);
+ split_pyramids(gr);
+ rearrange(gr);
+ }
- if (conformity == 1){
+ if (conformity >= 1){
init_markings(gr);
build_vertex_to_tetrahedra(gr);
build_vertex_to_pyramids(gr);
@@ -4971,6 +4976,43 @@ void PostOp::execute(GRegion* gr,int level, int conformity){
modify_surfaces(gr);
}
+void PostOp::init_markings_hex(GRegion* gr){
+ unsigned int i;
+ MElement* element;
+ markings.clear();
+ for(i=0;i<gr->getNumMeshElements();i++){
+ element = gr->getMeshElement(i);
+ if(eight(element)){
+ markings.insert(std::pair<MElement*,bool>(element,false));
+ }
+ }
+}
+
+void PostOp::init_markings_pri(GRegion* gr){
+ unsigned int i;
+ MElement* element;
+ markings.clear();
+ for(i=0;i<gr->getNumMeshElements();i++){
+ element = gr->getMeshElement(i);
+ if(six(element)){
+ markings.insert(std::pair<MElement*,bool>(element,false));
+ }
+ }
+}
+
+void PostOp::init_markings_pyr(GRegion* gr){
+ unsigned int i;
+ MElement* element;
+ markings.clear();
+ for(i=0;i<gr->getNumMeshElements();i++){
+ element = gr->getMeshElement(i);
+ if(five(element)){
+ markings.insert(std::pair<MElement*,bool>(element,false));
+ }
+ }
+}
+
+
void PostOp::init_markings(GRegion* gr){
unsigned int i;
MElement* element;
@@ -4985,6 +5027,223 @@ void PostOp::init_markings(GRegion* gr){
}
}
+void PostOp::split_hexahedra(GRegion* gr){
+ std::vector<MElement*> hexahedra;
+ std::vector<MHexahedron*> opt;
+ std::map<MElement*,bool>::iterator it;
+
+ hexahedra.clear();
+
+ for(size_t i=0;i<gr->getNumMeshElements();i++){
+ MElement *element = gr->getMeshElement(i);
+ if(eight(element)){
+ hexahedra.push_back(element);
+ }
+ }
+
+ for(size_t i=0;i<hexahedra.size();i++){
+ MElement *element = hexahedra[i];
+ MVertex *a = element->getVertex(0);
+ MVertex *b = element->getVertex(1);
+ MVertex *c = element->getVertex(2);
+ MVertex *d = element->getVertex(3);
+ MVertex *e = element->getVertex(4);
+ MVertex *f = element->getVertex(5);
+ MVertex *g = element->getVertex(6);
+ MVertex *h = element->getVertex(7);
+
+ bool conform=true;
+ conform &= (nonConformDiag(b,a,d,c,gr) == 0);
+ if (conform) conform &= (nonConformDiag(e,f,g,h,gr) == 0);
+ if (conform) conform &= (nonConformDiag(a,b,f,e,gr) == 0);
+ if (conform) conform &= (nonConformDiag(b,c,g,f,gr) == 0);
+ if (conform) conform &= (nonConformDiag(c,d,h,g,gr) == 0);
+ if (conform) conform &= (nonConformDiag(d,a,e,h,gr) == 0);
+ if (!conform){
+ double x = (a->x()+b->x()+c->x()+d->x()+e->x()+f->x()+g->x()+h->x())/8.0;
+ double y = (a->y()+b->y()+c->y()+d->y()+e->y()+f->y()+g->y()+h->y())/8.0;
+ double z = (a->z()+b->z()+c->z()+d->z()+e->z()+f->z()+g->z()+h->z())/8.0;
+ MVertex *mid = new MVertex(x,y,z,gr);
+ gr->addMeshVertex(mid);
+
+ MPyramid *temp = new MPyramid(a,b,c,d,mid);
+ gr->addPyramid(temp);
+ temp = new MPyramid(h,g,f,e,mid);
+ gr->addPyramid(temp);
+ temp = new MPyramid(e,f,b,a,mid);
+ gr->addPyramid(temp);
+ temp = new MPyramid(f,g,c,b,mid);
+ gr->addPyramid(temp);
+ temp = new MPyramid(g,h,d,c,mid);
+ gr->addPyramid(temp);
+ temp = new MPyramid(h,e,a,d,mid);
+ gr->addPyramid(temp);
+ it = markings.find(element);
+ it->second = true;
+ }
+ }
+
+ opt.clear();
+ opt.resize(gr->hexahedra.size());
+ opt = gr->hexahedra;
+ gr->hexahedra.clear();
+
+ for(size_t i=0;i<opt.size();i++){
+ MElement *element = (MElement*)(opt[i]);
+ it = markings.find(element);
+ if(it->second==0){
+ gr->hexahedra.push_back(opt[i]);
+ }
+ }
+}
+
+void PostOp::split_prisms(GRegion* gr){
+ std::vector<MElement*> prisms;
+ std::vector<MPrism*> opt;
+ std::map<MElement*,bool>::iterator it;
+
+ prisms.clear();
+
+ for(size_t i=0;i<gr->getNumMeshElements();i++){
+ MElement *element = gr->getMeshElement(i);
+ if(six(element)){
+ prisms.push_back(element);
+ }
+ }
+
+ for(size_t i=0;i<prisms.size();i++){
+ MElement *element = prisms[i];
+ MVertex *a = element->getVertex(0);
+ MVertex *b = element->getVertex(1);
+ MVertex *c = element->getVertex(2);
+ MVertex *d = element->getVertex(3);
+ MVertex *e = element->getVertex(4);
+ MVertex *f = element->getVertex(5);
+
+ pyramids1(a,d,f,c,gr);
+ pyramids1(a,b,e,d,gr);
+ pyramids1(b,c,f,e,gr);
+
+
+ bool conform=true;
+ conform &= (nonConformDiag(a,d,f,c,gr) == 0);
+ if (conform) conform &= (nonConformDiag(a,b,e,d,gr) == 0);
+ if (conform) conform &= (nonConformDiag(b,c,f,e,gr) == 0);
+ if (!conform){
+ double x = (a->x()+b->x()+c->x()+d->x()+e->x()+f->x())/6.0;
+ double y = (a->y()+b->y()+c->y()+d->y()+e->y()+f->y())/6.0;
+ double z = (a->z()+b->z()+c->z()+d->z()+e->z()+f->z())/6.0;
+ MVertex *mid = new MVertex(x,y,z,gr);
+ gr->addMeshVertex(mid);
+
+ MPyramid *temp = new MPyramid(c,f,d,a,mid);
+ gr->addPyramid(temp);
+ temp = new MPyramid(d,e,b,a,mid);
+ gr->addPyramid(temp);
+ temp = new MPyramid(e,f,c,b,mid);
+ gr->addPyramid(temp);
+ MTetrahedron *temp2 = new MTetrahedron(d,f,e,mid);
+ gr->addTetrahedron(temp2);
+ temp2 = new MTetrahedron(a,b,c,mid);
+ gr->addTetrahedron(temp2);
+ it = markings.find(element);
+ it->second = true;
+ }
+ }
+
+ opt.clear();
+ opt.resize(gr->prisms.size());
+ opt = gr->prisms;
+ gr->prisms.clear();
+
+ for(size_t i=0;i<opt.size();i++){
+ MElement *element = (MElement*)(opt[i]);
+ it = markings.find(element);
+ if(it->second==0){
+ gr->prisms.push_back(opt[i]);
+ }
+ }
+}
+
+void PostOp::split_pyramids(GRegion* gr){
+ std::vector<MElement*> pyramids;
+ std::vector<MPyramid*> opt;
+ std::map<MElement*,bool>::iterator it;
+
+ pyramids.clear();
+
+ for(size_t i=0;i<gr->getNumMeshElements();i++){
+ MElement *element = gr->getMeshElement(i);
+ if(five(element)){
+ pyramids.push_back(element);
+ }
+ }
+
+ for(size_t i=0;i<pyramids.size();i++){
+ MElement *element = pyramids[i];
+ MVertex *a = element->getVertex(0);
+ MVertex *b = element->getVertex(1);
+ MVertex *c = element->getVertex(2);
+ MVertex *d = element->getVertex(3);
+ MVertex *apex = element->getVertex(4);
+
+ int nDiag = nonConformDiag(a,b,c,d,gr);
+ if (nDiag == 1){
+ MTetrahedron *temp = new MTetrahedron(c,b,a,apex);
+ gr->addTetrahedron(temp);
+ temp = new MTetrahedron(c,a,d,apex);
+ gr->addTetrahedron(temp);
+ it = markings.find(element);
+ it->second = 1;
+ }else if (nDiag == 2){
+ MTetrahedron *temp = new MTetrahedron(b,a,d,apex);
+ gr->addTetrahedron(temp);
+ temp = new MTetrahedron(b,d,c,apex);
+ gr->addTetrahedron(temp);
+ it = markings.find(element);
+ it->second = 1;
+ }
+ }
+
+ opt.clear();
+ opt.resize(gr->pyramids.size());
+ opt = gr->pyramids;
+ gr->pyramids.clear();
+
+ for(size_t i=0;i<opt.size();i++){
+ MElement *element = (MElement*)(opt[i]);
+ it = markings.find(element);
+ if(it->second==0){
+ gr->pyramids.push_back(opt[i]);
+ }
+ }
+}
+
+
+int PostOp::nonConformDiag(MVertex* a,MVertex* b,MVertex* c,MVertex* d,GRegion* gr){
+ std::set<MElement*> diag1a;
+ std::set<MElement*> diag1b;
+ std::set<MElement*> diag2a;
+ std::set<MElement*> diag2b;
+
+ find_tetrahedra(a,b,c,diag1a);
+ find_tetrahedra(a,c,d,diag1b);
+ find_tetrahedra(b,c,d,diag2a);
+ find_tetrahedra(a,b,d,diag2b);
+ find_pyramids_from_tri(a,b,c,diag1a);
+ find_pyramids_from_tri(a,c,d,diag1b);
+ find_pyramids_from_tri(b,c,d,diag2a);
+ find_pyramids_from_tri(a,b,d,diag2b);
+ if(diag1a.size()==1 || diag1b.size()==1){
+ return 1;
+ }else if(diag2a.size()==1 || diag2b.size()==1){
+ return 2;
+ }
+ return 0;
+}
+
+
+
void PostOp::pyramids1(GRegion* gr){
unsigned int i;
MVertex *a,*b,*c,*d;
@@ -5199,8 +5458,8 @@ void PostOp::pyramids1(MVertex* a,MVertex* b,MVertex* c,MVertex* d,GRegion* gr){
gr->addPyramid(pyr);
it1->second = 1;
it2->second = 1;
- }else nbFAILVALIDITY++;
- }else nbFAILVERTEX++;
+ }
+ }
}
}
}
@@ -5421,9 +5680,6 @@ void PostOp::pyramids2(MVertex* a,MVertex* b,MVertex* c,MVertex* d,GRegion* gr,
MPyramid *pyr = new MPyramid(a,b,c,d,otherV[0]);
- if (otherV[0] != otherV[1])nbFAILVERTEX++;
- else if (!valid(pyr)) nbFAILVALIDITY++;
-
if (otherV[0] == otherV[1] && valid(pyr)){
estimate1 = estimate1 + tetrahedra.size() + 2*pyramids.size();
estimate2 = estimate2 + 1;
@@ -5432,7 +5688,7 @@ void PostOp::pyramids2(MVertex* a,MVertex* b,MVertex* c,MVertex* d,GRegion* gr,
//y = (diagA->y() + diagB->y() + otherV[0]->y())/3;
//z = (diagA->z() + diagB->z() + otherV[0]->z())/3;
- //W create a flat pyramid and let the optimizer fix it
+ //We create a flat pyramid and let the optimizer fix it
x = (diagA->x() + diagB->x())/2;
y = (diagA->y() + diagB->y())/2;
z = (diagA->z() + diagB->z())/2;
diff --git a/Mesh/yamakawa.h b/Mesh/yamakawa.h
index e4aa8c9e71bdb4d2369c0923d4e02f1734c22232..eadd202044b86915ff6924603885afbc1a16ed76 100644
--- a/Mesh/yamakawa.h
+++ b/Mesh/yamakawa.h
@@ -620,15 +620,13 @@ private:
std::multiset<Tuple> tuples;
std::set<MElement*> triangles;
- int nbFAILVERTEX;
- int nbFAILVALIDITY;
public:
PostOp();
~PostOp();
void execute(int, int);
- //level: 0: hex, 1: hex+prisms, 2: hex+prism+1-step pyramids, 3: hex+prism+2-steps pyramids
- //conformity= 0=nonconforming, 1=conformity using trihedra
+ //level - 0: hex, 1: hex+prisms, 2: hex+prism+pyramids
+ //conformity - 0: nonconforming, 1: trihedra, 2: pyramids+trihedra, 3:pyramids+hexPrismSplit+trihedra, 4:hexPrismSplit+trihedra
void execute(GRegion*,int level, int conformity);
inline int get_nb_hexahedra()const{return nbr8;};
@@ -637,9 +635,16 @@ public:
inline double get_vol_elements()const{return vol;};
void init_markings(GRegion*);
+ void init_markings_hex(GRegion*);
+ void init_markings_pri(GRegion*);
+ void init_markings_pyr(GRegion*);
void pyramids1(GRegion*);
void pyramids2(GRegion*, bool allowNonConforming=false);
void trihedra(GRegion*);
+ void split_hexahedra(GRegion*);
+ void split_prisms(GRegion*);
+ void split_pyramids(GRegion*);
+ int nonConformDiag(MVertex* a,MVertex* b,MVertex* c,MVertex* d,GRegion* gr);
void pyramids1(MVertex*,MVertex*,MVertex*,MVertex*,GRegion*);
void pyramids2(MVertex*,MVertex*,MVertex*,MVertex*,GRegion*, bool allowNonConforming);
void trihedra(MVertex*,MVertex*,MVertex*,MVertex*,GRegion*);