diff --git a/Mesh/CMakeLists.txt b/Mesh/CMakeLists.txt
index d6db42e7055d8d3ea21d7a06c2b53e320b7c0072..0a27aec8a8bcccc21f359d84159ebb2cd7e29333 100644
--- a/Mesh/CMakeLists.txt
+++ b/Mesh/CMakeLists.txt
@@ -20,6 +20,7 @@ set(SRC
meshGRegionDelaunayInsertion.cpp meshGRegionTransfinite.cpp
meshGRegionExtruded.cpp meshGRegionCarveHole.cpp
meshGRegionLocalMeshMod.cpp meshGRegionMMG3D.cpp
+meshGRegionRelocateVertex.cpp
meshMetric.cpp
BackgroundMeshTools.cpp
BackgroundMesh.cpp
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 430a159c6970cdd4df2e1f354b6603f3256befc7..20e13256c1bcd8ee10d21c50b9d8ef9c09aa8c18 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -13,6 +13,7 @@
#include "boundaryLayersData.h"
#include "meshGRegionBoundaryRecovery.h"
#include "meshGRegionDelaunayInsertion.h"
+#include "meshGRegionRelocateVertex.h"
#include "GModel.h"
#include "GRegion.h"
#include "GFace.h"
@@ -120,112 +121,6 @@ class splitQuadRecovery {
}
};
-void printVoronoi(GRegion *gr, std::set<SPoint3> &candidates)
-{
- std::vector<MTetrahedron*> elements = gr->tetrahedra;
- std::list<GFace*> allFaces = gr->faces();
-
- //building maps
- std::map<MVertex*, std::set<MTetrahedron*> > node2Tet;
- std::map<MFace, std::vector<MTetrahedron*> , Less_Face> face2Tet;
- for(unsigned int i = 0; i < elements.size(); i++){
- MTetrahedron *ele = elements[i];
- for (int j=0; j<4; j++){
- MVertex *v = ele->getVertex(j);
- std::map<MVertex*, std::set<MTetrahedron*> >::iterator itmap = node2Tet.find(v);
- if (itmap == node2Tet.end()){
- std::set<MTetrahedron*> oneTet;
- oneTet.insert(ele);
- node2Tet.insert(std::make_pair(v, oneTet));
- }
- else{
- std::set<MTetrahedron*> allTets = itmap->second;
- allTets.insert(allTets.begin(), ele);
- itmap->second = allTets;
- }
- }
- for (int j = 0; j < 4; j++){
- MFace f = ele->getFace(j);
- std::map<MFace, std::vector<MTetrahedron*>, Less_Face >::iterator itmap =
- face2Tet.find(f);
- if (itmap == face2Tet.end()){
- std::vector<MTetrahedron*> oneTet;
- oneTet.push_back(ele);
- face2Tet.insert(std::make_pair(f, oneTet));
- }
- else{
- std::vector<MTetrahedron*> allTets = itmap->second;
- allTets.insert(allTets.begin(), ele);
- itmap->second = allTets;
- }
- }
- }
-
- //print voronoi poles
- //smooth_normals *snorm = gr->model()->normals;
- FILE *outfile = Fopen("nodes.pos", "w");
- if(outfile) fprintf(outfile, "View \"Voronoi poles\" {\n");
- std::map<MVertex*, std::set<MTetrahedron*> >::iterator itmap = node2Tet.begin();
- for(; itmap != node2Tet.end(); itmap++){
- //MVertex *v = itmap->first;
- std::set<MTetrahedron*> allTets = itmap->second;
- std::set<MTetrahedron*>::const_iterator it = allTets.begin();
- MTetrahedron *poleTet = *it;
- double maxRadius = poleTet->getCircumRadius();
- for(; it != allTets.end(); it++){
- double radius = (*it)->getCircumRadius();
- if (radius > maxRadius){
- maxRadius = radius;
- poleTet = *it;
- }
- }
- //if (v->onWhat()->dim() == 2 ){
- SPoint3 pc = poleTet->circumcenter();
- //double nx,ny,nz;
- // SVector3 vN = snorm->get(v->x(), v->y(), v->z(), nx,ny,nz);
- // SVector3 pcv(pc.x()-nx, pc.y()-ny,pc.z()-nz);
- // printf("nx=%g ny=%g nz=%g dot=%g \n", nx,ny,nz, dot(vN, pcv));
- // if ( dot(vN, pcv) > 0.0 )
- double x[3] = {pc.x(), pc.y(), pc.z()};
- double uvw[3];
- poleTet->xyz2uvw(x, uvw);
- //bool inside = poleTet->isInside(uvw[0], uvw[1], uvw[2]);
- //if (inside){
- if(outfile)
- fprintf(outfile,"SP(%g,%g,%g) {%g};\n", pc.x(), pc.y(), pc.z(), maxRadius);
- candidates.insert(pc);
- //}
- //}
- }
- if(outfile){
- fprintf(outfile,"};\n");
- fclose(outfile);
- }
-
- // print scalar lines
- FILE *outfile2 = Fopen("edges.pos", "w");
- if(outfile2) fprintf(outfile2, "View \"Voronoi edges\" {\n");
- std::map<MFace, std::vector<MTetrahedron*> , Less_Face >::iterator itmap2 = face2Tet.begin();
- for(; itmap2 != face2Tet.end(); itmap2++){
- std::vector<MTetrahedron*> allTets = itmap2->second;
- if (allTets.size() != 2 ) continue;
- SPoint3 pc1 = allTets[0]->circumcenter();
- SPoint3 pc2 = allTets[1]->circumcenter();
- //std::set<SPoint3>::const_iterator it1 = candidates.find(pc1);
- //std::set<SPoint3>::const_iterator it2 = candidates.find(pc2);
- //if( it1 != candidates.end() || it2 != candidates.end())
- if(outfile2)
- fprintf(outfile2,"SL(%g,%g,%g,%g,%g,%g) {%g,%g};\n",
- pc1.x(), pc1.y(), pc1.z(), pc2.x(), pc2.y(), pc2.z(),
- allTets[0]->getCircumRadius(),allTets[1]->getCircumRadius());
- }
- if(outfile2){
- fprintf(outfile2,"};\n");
- fclose(outfile2);
- }
-}
-
-
void getBoundingInfoAndSplitQuads(GRegion *gr,
std::map<MFace,GEntity*,Less_Face> &allBoundingFaces,
std::set<MVertex*> &allBoundingVertices,
@@ -550,613 +445,6 @@ static void addOrRemove(const MFace &f,
}
-/*
- here, we have a list of elements that actually do not form a mesh
- --> a set boundary layer elements (prism, hexes, pyramids (and soon tets)
- --> a set of tetrahedra that respect the external boundary of the BL mesh,
- yet possiblty containing elements that are on the other side of the boundary
- and therefore overlapping those elements. We have to extract the good ones !
-
-
-
-*/
-
-static bool AssociateElementsToModelRegionWithBoundaryLayers (GRegion *gr,
- std::vector<MTetrahedron*> &tets,
- std::vector<MHexahedron*> &hexes,
- std::vector<MPrism*> &prisms,
- std::vector<MPyramid*> &pyramids,
- splitQuadRecovery & sqr)
-{
- std::set<MElement*> all;
- all.insert(hexes.begin(),hexes.end());
- all.insert(prisms.begin(),prisms.end());
- all.insert(pyramids.begin(),pyramids.end());
- // start with one BL element !
- MElement *FIRST = all.size() ? *(all.begin()) : 0;
- if (!FIRST) return true;
- all.insert(tets.begin(),tets.end());
-
- // printf("coucou1 %d eleemnts\n",all.size());
- fs_cont search;
- buildFaceSearchStructure(gr->model(), search);
-
- // create the graph face 2 elements for the region
- std::map<MFace,std::pair<MElement*,MElement*> ,Less_Face> myGraph;
-
- for (std::set<MElement*>::iterator it = all.begin(); it != all.end(); ++it){
- MElement *t = *it;
- const int nbf = t->getNumFaces();
- for (int j=0;j<nbf;j++){
- MFace f = t->getFace(j);
- std::map<MFace,std::pair<MElement*,MElement*>,Less_Face>::iterator itf = myGraph.find(f);
- if (itf == myGraph.end())myGraph.insert (std::make_pair (f, std::make_pair (t,(MElement*)0)));
- else {
- // what to do ??????
- // two tets and one prism --> the prism should be
- // geometrically on the other side of the tet
- if (itf->second.second) {
- MElement *prism=0, *t1=0, *t2=0;
- if (itf->second.second->getType () == TYPE_PRI || itf->second.second->getType () == TYPE_PYR) {
- prism = itf->second.second;
- t1 = itf->second.first;
- t2 = t;
- }
- else if (itf->second.first->getType () == TYPE_PRI || itf->second.first->getType () == TYPE_PYR) {
- prism = itf->second.first;
- t1 = itf->second.second;
- t2 = t;
- }
- else if (t->getType () == TYPE_PRI || t->getType () == TYPE_PYR) {
- prism = t;
- t1 = itf->second.second;
- t2 = itf->second.first;
- }
- else {
- printf("types %d %d %d\n",t->getType () ,itf->second.first->getType (),itf->second.second->getType () );
- }
- if (!t1 || !t2 || !prism){
- gr->model()->writeMSH("ooops.msh");
- Msg::Error ("Wrong BL configuration");
- return false;
- }
- SVector3 n = f.normal();
- SPoint3 c = f.barycenter();
- MVertex *v_prism = 0, *v_t1 = 0, *v_t2 = 0;
- for (int i=0;i<4;i++){
- if (t1->getVertex(i) != f.getVertex(0) &&
- t1->getVertex(i) != f.getVertex(1) &&
- t1->getVertex(i) != f.getVertex(2))v_t1 = t1->getVertex(i);
- if (t2->getVertex(i) != f.getVertex(0) &&
- t2->getVertex(i) != f.getVertex(1) &&
- t2->getVertex(i) != f.getVertex(2))v_t2 = t2->getVertex(i);
- }
- for (int i=0;i<prism->getNumVertices();i++){
- if (prism->getVertex(i) != f.getVertex(0) &&
- prism->getVertex(i) != f.getVertex(1) &&
- prism->getVertex(i) != f.getVertex(2)) v_prism = prism->getVertex(i);
- }
- SVector3 vf1 (v_t1->x() - c.x(),v_t1->y() - c.y(),v_t1->z() - c.z());
- SVector3 vf2 (v_t2->x() - c.x(),v_t2->y() - c.y(),v_t2->z() - c.z());
- SVector3 vfp (v_prism->x() - c.x(),v_prism->y() - c.y(),v_prism->z() - c.z());
- // printf("%12.5E %12.5E%12.5E\n",dot(vf1,n),dot(vf2,n),dot(vfp,n));
- if (dot (vf1,n) * dot (vfp,n) > 0){itf->second.first = prism;itf->second.second=t2; /*delete t1;*/}
- else if (dot (vf2,n) * dot (vfp,n) > 0){itf->second.first = prism;itf->second.second=t1; /*delete t2;*/}
- // else if (dot (vf2,vfp) > 0){itf->second.first = prism;itf->second.second=t2;}
- else Msg::Fatal("Impossible Geom Config");
- }
- else itf->second.second = t;
- }
- }
- }
-
- std::stack<MElement*> myStack;
- std::set<MElement*> connected;
- std::set<GFace*> faces_bound;
- myStack.push(FIRST);
- while (!myStack.empty()){
- FIRST = myStack.top();
- myStack.pop();
- connected.insert(FIRST);
- for (int i=0;i<FIRST->getNumFaces();i++){
- MFace f = FIRST->getFace(i);
- std::map<MFace, MVertex*, Less_Face>::iterator it = sqr._invmap.find(f);
- GFace* gfound = 0;
- if (it != sqr._invmap.end()){
- gfound = (GFace*)it->second->onWhat();
- // one pyramid is useless because one element with a quad face impacts the
- // boundary of the domain.
- sqr._toDelete.insert(f);
- }
- else gfound = findInFaceSearchStructure (f,search);
- if (!gfound){
- std::map<MFace,std::pair<MElement*,MElement*>,Less_Face>::iterator
- itf = myGraph.find(f);
- MElement *t_neigh = itf->second.first == FIRST ?
- itf->second.second : itf->second.first;
- if (!t_neigh)printf("oulalalalalalalala %d vertices\n",f.getNumVertices());
- if (connected.find(t_neigh) == connected.end())myStack.push(t_neigh);
- }
- else {
- // if (faces_bound.find(gfound) == faces_bound.end())printf("face %d\n",gfound->tag());
- faces_bound.insert(gfound);
- }
- }
- }
-
- // printf ("found a set of %d elements that are connected with %d bounding faces\n",connected.size(),faces_bound.size());
- GRegion *myGRegion = getRegionFromBoundingFaces(gr->model(), faces_bound);
- // printf("REGION %d %d\n",myGRegion->tag(),gr->tag());
- if (myGRegion == gr){
- // printf("one region %d is found : %d elements\n",myGRegion->tag(),connected.size());
- myGRegion->tetrahedra.clear();
- for (std::set<MElement*>::iterator it = connected.begin(); it != connected.end(); ++it){
- MElement *t = *it;
- std::set<MVertex*> _mesh_vertices;
- for (int i=0;i<t->getNumVertices();i++){
- MVertex *_v = t->getVertex(i);
- if (_v->onWhat() == gr){
- _mesh_vertices.insert(_v);
- }
- }
- // myGRegion->mesh_vertices.insert(myGRegion->mesh_vertices.end(),_mesh_vertices.begin(),_mesh_vertices.end());
-
- if (t->getType() == TYPE_TET)
- myGRegion->tetrahedra.push_back((MTetrahedron*)t);
- else if (t->getType() == TYPE_HEX)
- myGRegion->hexahedra.push_back((MHexahedron*)t);
- else if (t->getType() == TYPE_PRI)
- myGRegion->prisms.push_back((MPrism*)t);
- else if (t->getType() == TYPE_PYR)
- myGRegion->pyramids.push_back((MPyramid*)t);
- }
- }
- else {
- return false;
- }
- return true;
-}
-
-static int getWedge(BoundaryLayerColumns* _columns, MVertex *v1, MVertex *v2,
- int indicesVert1 [], int indicesVert2 [])
-{
- int N1 = _columns->getNbColumns(v1) ;
- int N2 = _columns->getNbColumns(v2) ;
- int fanSize = 4;
- int NW1 = 0;
- int NW2 = 0;
- for (int i=0;i<N1;i++){
- const BoundaryLayerData & c1 = _columns->getColumn(v1,i);
- if (c1._joint.size())NW1++;
- }
- for (int i=0;i<N2;i++){
- const BoundaryLayerData & c2 = _columns->getColumn(v2,i);
- if (c2._joint.size())NW2++;
- }
-
- std::map<int,int> one2two;
- for (int i=0;i<NW1;i++){
- const BoundaryLayerData & c1 = _columns->getColumn(v1,i);
- for (int j=0;j<NW2;j++){
- const BoundaryLayerData & c2 = _columns->getColumn(v2,j);
- for (unsigned int k=0;k<c2._joint.size();k++){
- if (std::find(c1._joint.begin(),c1._joint.end(),c2._joint[k]) !=
- c1._joint.end()) {
- one2two[i] = j;
- }
- }
- }
- }
-
- // printf("%d %d %d %d \n",N1,N2,NW1,NW2);
- // for (std::map<int,int>::iterator it = one2two.begin(); it != one2two.end(); it++){
- // printf("one2two[%d] = %d\n",it->first,it->second);
- // }
- if (one2two.size() != 2)return 0;
-
- int vert1Start,vert1End;
- int vert2Start,vert2End;
- std::map<int,int>::iterator it = one2two.begin();
- vert1Start = it->first;
- vert2Start = it->second;
- ++it;
- vert1End = it->first;
- vert2End = it->second;
-
-
- int INDEX1 = 0, count = 0;
- for (int i=0;i<NW1;i++){
- for (int j=i+1;j<NW1;j++){
- if ((vert1Start == i && vert1End == j) ||
- (vert1Start == j && vert1End == i))
- {
- INDEX1 = count;
- }
- count++;
- }
- }
- int INDEX2 = 0;
- count = 0;
- for (int i=0;i<NW2;i++){
- for (int j=i+1;j<NW2;j++){
- if ((vert2Start == i && vert2End == j) ||
- (vert2Start == j && vert2End == i))
- {
- INDEX2 = count;
- }
- count++;
- }
- }
-
- int indexVert1Start = NW1 + fanSize * ( 0 + INDEX1);
- int indexVert1End = NW1 + fanSize * ( 1 + INDEX1);
-
- int indexVert2Start = NW2 + fanSize * ( 0 + INDEX2);
- int indexVert2End = NW2 + fanSize * ( 1 + INDEX2);
-
- indicesVert1[0] = vert1Start;
- int k=1;
- for (int i=indexVert1Start;i< indexVert1End;++i)indicesVert1[k++] = i;
- indicesVert1[fanSize+1] = vert1End;
-
- indicesVert2[0] = vert2Start;
- k = 1;
- if (indexVert2End > indexVert2Start){
- for (int i=indexVert2Start;i< indexVert2End;++i)indicesVert2[k++] = i;
- }
- else {
- for (int i=indexVert2Start-1;i<= indexVert2End;--i)indicesVert2[k++] = i;
- }
- indicesVert2[fanSize+1] = vert2End;
-
-
- // printf("%d %d %d %d %d %d %d %d\n",vert1Start,vert1End,vert2Start,vert2End,indexVert1Start,indexVert1End,indexVert2Start,indexVert2End);
- // return 0;
-
- return fanSize + 2;
-}
-
-
-static bool modifyInitialMeshForTakingIntoAccountBoundaryLayers(GRegion *gr, splitQuadRecovery & sqr)
-{
- if (getBLField(gr->model())) insertVerticesInRegion(gr,-1);
- if (!buildAdditionalPoints3D (gr)) return false;
- BoundaryLayerColumns* _columns = gr->getColumns();
- std::map<MFace,MElement*,Less_Face> bfaces;
-
- std::vector<MPrism*> blPrisms;
- std::vector<MHexahedron*> blHexes;
- std::vector<MPyramid*> blPyrs;
- std::list<GFace*> faces = gr->faces();
-
- std::list<GFace*> embedded_faces = gr->embeddedFaces();
- faces.insert(faces.begin(), embedded_faces.begin(),embedded_faces.end());
- std::set<MVertex*> verts;
-
- {
- std::list<GFace*>::iterator itf = faces.begin();
- while(itf != faces.end()){
- for(unsigned int i = 0; i< (*itf)->getNumMeshElements(); i++){
- MElement *e = (*itf)->getMeshElement(i);
- addOrRemove (e->getFace(0),0,bfaces,sqr);
- }
- ++itf;
- }
- }
-
- std::list<GFace*>::iterator itf = faces.begin();
- while(itf != faces.end()){
- for(unsigned int i = 0; i< (*itf)->triangles.size(); i++){
- MVertex *v1 = (*itf)->triangles[i]->getVertex(0);
- MVertex *v2 = (*itf)->triangles[i]->getVertex(1);
- MVertex *v3 = (*itf)->triangles[i]->getVertex(2);
- MFace dv (v1,v2,v3);
- for (unsigned int SIDE = 0 ; SIDE < _columns->_normals3D.count(dv); SIDE ++){
- faceColumn fc = _columns->getColumns(*itf,v1, v2, v3, SIDE);
- const BoundaryLayerData & c1 = fc._c1;
- const BoundaryLayerData & c2 = fc._c2;
- const BoundaryLayerData & c3 = fc._c3;
- int N = std::min(c1._column.size(),std::min(c2._column.size(),c3._column.size()));
-
- // double distMax = getDistMax (v1, v2, v3, c1._n, c2._n, c3._n);
- MFace f_low (v1,v2,v3);
- SVector3 n_low = f_low.normal();
- // printf("%d Layers\n",N);
- std::vector<MElement*> myCol;
- for (int l=0;l < N ;++l){
- MVertex *v11,*v12,*v13,*v21,*v22,*v23;
- v21 = c1._column[l];
- v22 = c2._column[l];
- v23 = c3._column[l];
- if (l == 0){
- v11 = v1;
- v12 = v2;
- v13 = v3;
- }
- else {
- v11 = c1._column[l-1];
- v12 = c2._column[l-1];
- v13 = c3._column[l-1];
- }
- MFace f_up (v21,v22,v23);
- SVector3 n_up = f_up.normal();
- double dotProd = dot(n_up,n_low);
- MPrism *prism = new MPrism(v11,v12,v13,v21,v22,v23);
- if (dotProd > 0.2 && prism->skewness() > 0.1){
- blPrisms.push_back(prism);
- myCol.push_back(prism);
- }
- else {
- delete prism;
- l = N+1;
- }
- }
- if (!myCol.empty()){
- for (unsigned int l=0;l<myCol.size();l++)_columns->_toFirst[myCol[l]] = myCol[0];
- _columns->_elemColumns[myCol[0]] = myCol;
- }
- }
- }
- ++itf;
- }
-
- std::set<MEdge,Less_Edge> edges;
- {
- std::list<GEdge*> gedges = gr->edges();
- for (std::list<GEdge*>::iterator it = gedges.begin(); it != gedges.end() ; ++it){
- for (unsigned int i=0;i<(*it)->lines.size();++i){
- edges.insert(MEdge((*it)->lines[i]->getVertex(0),(*it)->lines[i]->getVertex(1)));
- }
- }
- }
-
- // now treat the Wedges
- // we have to know the two target GFaces that are concerned with a GEdge
- std::set<MEdge,Less_Edge>::iterator ite = edges.begin();
- while(ite != edges.end()){
- MEdge e = *ite;
- MVertex *v1 = e.getVertex(0);
- MVertex *v2 = e.getVertex(1);
- if (v1 != v2){
- int indices1[256];
- int indices2[256];
- int NbW = getWedge (_columns, v1, v2, indices1,indices2);
- for (int i=0;i<NbW-1;i++){
- int i11 = indices1[i];
- int i12 = indices1[i+1];
- int i21 = indices2[i];
- int i22 = indices2[i+1];
- BoundaryLayerData c11 = _columns->getColumn(v1,i11);
- BoundaryLayerData c12 = _columns->getColumn(v1,i12);
- BoundaryLayerData c21 = _columns->getColumn(v2,i21);
- BoundaryLayerData c22 = _columns->getColumn(v2,i22);
- int N = std::min(c11._column.size(),
- std::min(c12._column.size(),
- std::min(c21._column.size(), c22._column.size())));
- std::vector<MElement*> myCol;
- for (int l=0;l < N ;++l){
- MVertex *v11,*v12,*v13,*v14;
- MVertex *v21,*v22,*v23,*v24;
- v21 = c11._column[l];
- v22 = c12._column[l];
- v23 = c22._column[l];
- v24 = c21._column[l];
- if (l == 0){
- v11 = v12 = v1;
- v13 = v14 = v2;
- }
- else {
- v11 = c11._column[l-1];
- v12 = c12._column[l-1];
- v13 = c22._column[l-1];
- v14 = c21._column[l-1];
- }
-
- if (l == 0){
- MPrism *prism = new MPrism(v12,v21,v22,v13,v24,v23);
- // store the layer the element belongs
- myCol.push_back(prism);
-
- blPrisms.push_back(prism);
- }
- else {
- MHexahedron *hex = new MHexahedron(v11,v12,v13,v14,v21,v22,v23,v24);
- // store the layer the element belongs
- myCol.push_back(hex);
- blHexes.push_back(hex);
- }
- }
- if (!myCol.empty()){
- for (unsigned int l=0;l<myCol.size();l++)_columns->_toFirst[myCol[l]] = myCol[0];
- _columns->_elemColumns[myCol[0]] = myCol;
- }
- }
- }
- ++ite;
- }
- // ------------------------------------------------------------------------------------
- // FIXME : NOT 100 % CORRECT
- // filterOverlappingElements (blPrisms,blHexes,_columns->_elemColumns,_columns->_toFirst);
- // ------------------------------------------------------------------------------------
- {
- FILE *ff2 = fopen ("tato3D.pos","w");
- if(ff2){
- fprintf(ff2,"View \" \"{\n");
- for (unsigned int i = 0; i < blPrisms.size();i++){
- blPrisms[i]->writePOS(ff2,1,0,0,0,0,0);
- }
- for (unsigned int i = 0; i < blHexes.size();i++){
- blHexes[i]->writePOS(ff2,1,0,0,0,0,0);
- }
- fprintf(ff2,"};\n");
- fclose(ff2);
- }
- }
-
- for (unsigned int i = 0; i < blPrisms.size();i++){
- for (unsigned int j=0;j<5;j++)
- addOrRemove(blPrisms[i]->getFace(j),blPrisms[i],bfaces,sqr);
- for (int j = 0; j < 6; j++)
- if(blPrisms[i]->getVertex(j)->onWhat() == gr)verts.insert(blPrisms[i]->getVertex(j));
- }
- for (unsigned int i = 0; i < blHexes.size();i++){
- for (unsigned int j=0;j<6;j++)
- addOrRemove(blHexes[i]->getFace(j),blHexes[i],bfaces, sqr);
- for (int j = 0; j < 8; j++)
- if(blHexes[i]->getVertex(j)->onWhat() == gr)verts.insert(blHexes[i]->getVertex(j));
- }
-
- discreteFace *nf = new discreteFace(gr->model(), 444444);
- gr->model()->add (nf);
- std::list<GFace*> hop;
- std::map<MFace,MElement*,Less_Face>::iterator it = bfaces.begin();
-
- FILE *ff = fopen ("toto3D.pos","w");
- if(ff) fprintf(ff,"View \" \"{\n");
- for (; it != bfaces.end(); ++it){
- if (it->first.getNumVertices() == 3){
- nf->triangles.push_back(new MTriangle (it->first.getVertex(0),it->first.getVertex(1),it->first.getVertex(2)));
- if(ff)
- fprintf(ff,"ST (%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1};\n",
- it->first.getVertex(0)->x(),it->first.getVertex(0)->y(),it->first.getVertex(0)->z(),
- it->first.getVertex(1)->x(),it->first.getVertex(1)->y(),it->first.getVertex(1)->z(),
- it->first.getVertex(2)->x(),it->first.getVertex(2)->y(),it->first.getVertex(2)->z());
- }
- else {
-
- // we have a quad face --> create a pyramid
-
- MElement *e = it->second;
-
- // compute the center of the face;
- SPoint3 center (0.25*(it->first.getVertex(0)->x()+it->first.getVertex(1)->x()+it->first.getVertex(2)->x()+it->first.getVertex(3)->x()),
- 0.25*(it->first.getVertex(0)->y()+it->first.getVertex(1)->y()+it->first.getVertex(2)->y()+it->first.getVertex(3)->y()),
- 0.25*(it->first.getVertex(0)->z()+it->first.getVertex(1)->z()+it->first.getVertex(2)->z()+it->first.getVertex(3)->z()));
- // compute the center of the opposite face;
- SPoint3 opposite (0,0,0);
- int counter=0;
- for (int i=0;i<e->getNumVertices();i++){
- MVertex *vv = e->getVertex(i);
- bool found = false;
- for (int j=0;j<4;j++){
- MVertex *ww = it->first.getVertex(j);
- if (ww == vv)found = true;
- }
- if (!found){
- counter ++;
- opposite += SPoint3(vv->x(),vv->y(),vv->z());
- }
- }
- // printf("counter = %d\n",counter);
- if(counter)
- opposite /= (double)counter;
-
- SVector3 dir = center - opposite;
- MTriangle temp (it->first.getVertex(0),it->first.getVertex(1),it->first.getVertex(2));
- double D = temp.minEdge();
- dir.normalize();
- const double eps = D*1.e-2;
- MVertex *newv = new MVertex (center.x() + eps * dir.x(),center.y() + eps * dir.y(), center.z() + eps * dir.z(), gr);
-
- MPyramid *pyr = new MPyramid (it->first.getVertex(0),it->first.getVertex(1),it->first.getVertex(2),it->first.getVertex(3),newv);
- verts.insert(newv);
- blPyrs.push_back(pyr);
-
- nf->triangles.push_back(new MTriangle (it->first.getVertex(0),it->first.getVertex(1),newv));
- nf->triangles.push_back(new MTriangle (it->first.getVertex(1),it->first.getVertex(2),newv));
- nf->triangles.push_back(new MTriangle (it->first.getVertex(2),it->first.getVertex(3),newv));
- nf->triangles.push_back(new MTriangle (it->first.getVertex(3),it->first.getVertex(0),newv));
- if(ff){
- fprintf(ff,"ST (%g,%g,%g,%g,%g,%g,%g,%g,%g){2,2,2};\n",
- it->first.getVertex(0)->x(),it->first.getVertex(0)->y(),it->first.getVertex(0)->z(),
- it->first.getVertex(1)->x(),it->first.getVertex(1)->y(),it->first.getVertex(1)->z(),
- newv->x(),newv->y(),newv->z());
-
- fprintf(ff,"ST (%g,%g,%g,%g,%g,%g,%g,%g,%g){2,2,2};\n",
- it->first.getVertex(1)->x(),it->first.getVertex(1)->y(),it->first.getVertex(1)->z(),
- it->first.getVertex(2)->x(),it->first.getVertex(2)->y(),it->first.getVertex(2)->z(),
- newv->x(),newv->y(),newv->z());
-
- fprintf(ff,"ST (%g,%g,%g,%g,%g,%g,%g,%g,%g){2,2,2};\n",
- it->first.getVertex(2)->x(),it->first.getVertex(2)->y(),it->first.getVertex(2)->z(),
- it->first.getVertex(3)->x(),it->first.getVertex(3)->y(),it->first.getVertex(3)->z(),
- newv->x(),newv->y(),newv->z());
- fprintf(ff,"ST (%g,%g,%g,%g,%g,%g,%g,%g,%g){2,2,2};\n",
- it->first.getVertex(3)->x(),it->first.getVertex(3)->y(),it->first.getVertex(3)->z(),
- it->first.getVertex(0)->x(),it->first.getVertex(0)->y(),it->first.getVertex(0)->z(),
- newv->x(),newv->y(),newv->z());
-
- fprintf(ff,"SQ (%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){3,3,3,3};\n",
- it->first.getVertex(0)->x(),it->first.getVertex(0)->y(),it->first.getVertex(0)->z(),
- it->first.getVertex(1)->x(),it->first.getVertex(1)->y(),it->first.getVertex(1)->z(),
- it->first.getVertex(2)->x(),it->first.getVertex(2)->y(),it->first.getVertex(2)->z(),
- it->first.getVertex(3)->x(),it->first.getVertex(3)->y(),it->first.getVertex(3)->z());
- }
- }
- }
- if(ff){
- fprintf(ff,"};\n");
- fclose(ff);
- }
-
- printf("discrete face with %d %d elems\n", (int)nf->triangles.size(),
- (int)nf->quadrangles.size());
- hop.push_back(nf);
-
- for(unsigned int i = 0; i < gr->tetrahedra.size(); i++) delete gr->tetrahedra[i];
- gr->tetrahedra.clear();
-
- gr->set(hop);
- CreateAnEmptyVolumeMesh(gr);
- printf("%d tets\n", (int)gr->tetrahedra.size());
- deMeshGFace _kill;
- _kill (nf);
- //<<<<<<< .mine
- gr->model()->remove(nf);
- // delete nf;
- //=======
- //>>>>>>> .r18259
- delete nf;
-
- gr->set(faces);
- gr->mesh_vertices.insert(gr->mesh_vertices.begin(),verts.begin(),verts.end());
-
- gr->model()->writeMSH("BL_start.msh");
-
- AssociateElementsToModelRegionWithBoundaryLayers (gr, gr->tetrahedra , blHexes, blPrisms, blPyrs, sqr);
-
- gr->model()->writeMSH("BL_start2.msh");
-
- return true;
-}
-
-void _relocateVertex(MVertex *ver,
- const std::vector<MElement*> <)
-{
- if(ver->onWhat()->dim() != 3) return;
- double x = 0, y=0, z=0;
- int N = 0;
- bool isPyramid = false;
- for(unsigned int i = 0; i < lt.size(); i++){
- double XCG=0,YCG=0,ZCG=0;
- for (int j=0;j<lt[i]->getNumVertices();j++){
- XCG += lt[i]->getVertex(j)->x();
- YCG += lt[i]->getVertex(j)->y();
- ZCG += lt[i]->getVertex(j)->z();
- }
- x += XCG;
- y += YCG;
- z += ZCG;
- if (lt[i]->getNumVertices() == 5) isPyramid = true;
- N += lt[i]->getNumVertices();
- }
- if (isPyramid){
- ver->x() = x / N;
- ver->y() = y / N;
- ver->z() = z / N;
- }
-}
-
#if defined(HAVE_TETGEN)
bool CreateAnEmptyVolumeMesh(GRegion *gr)
{
@@ -1282,8 +570,6 @@ void MeshDelaunayVolumeTetgen(std::vector<GRegion*> ®ions)
// restore the initial set of faces
gr->set(faces);
- bool _BL = modifyInitialMeshForTakingIntoAccountBoundaryLayers(gr,sqr);
-
// now do insertion of points
if(CTX::instance()->mesh.algo3d == ALGO_3D_FRONTAL_DEL)
bowyerWatsonFrontalLayers(gr, false);
@@ -1295,40 +581,12 @@ void MeshDelaunayVolumeTetgen(std::vector<GRegion*> ®ions)
else{
int nbvertices_filler = (old_algo_hexa()) ? Filler::get_nbr_new_vertices() : Filler3D::get_nbr_new_vertices();
if(!nbvertices_filler && !LpSmoother::get_nbr_interior_vertices()){
- insertVerticesInRegion(gr,2000000000,!_BL);
+ insertVerticesInRegion(gr,2000000000,true);
}
}
- //emi test frame field
- // int NumSmooth = 10;//CTX::instance()->mesh.smoothCrossField
- // std::cout << "NumSmooth = " << NumSmooth << std::endl;
- // if(NumSmooth && (gr->dim() == 3)){
- // double scale = gr->bounds().diag()*1e-2;
- // Frame_field::initRegion(gr,NumSmooth);
- // Frame_field::saveCrossField("cross0.pos",scale);
- // //Frame_field::smoothRegion(gr,NumSmooth);
- // //Frame_field::saveCrossField("cross1.pos",scale);
- // GFace *gf = GModel::current()->getFaceByTag(2);
- // Frame_field::continuousCrossField(gr,gf);
- // Frame_field::saveCrossField("cross2.pos",scale);
- // }
- // Frame_field::init_region(gr);
- // Frame_field::clear();
- // exit(1);
- //fin test emi
-
if (sqr.buildPyramids (gr->model())){
- // relocate vertices if pyramids
- for(unsigned int k = 0; k < regions.size(); k++){
- v2t_cont adj;
- buildVertexToElement(regions[k]->tetrahedra, adj);
- buildVertexToElement(regions[k]->pyramids, adj);
- v2t_cont::iterator it = adj.begin();
- while (it != adj.end()){
- _relocateVertex( it->first, it->second);
- ++it;
- }
- }
+ RelocateVertices (regions);
}
#endif
}
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index 4d252f5c96716f285c481a02c41c977263fa0a33..adc94460b077a9a9e504c4c74f8de1556b3e8444 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -815,7 +815,6 @@ void adaptMeshGRegion::operator () (GRegion *gr)
}
}
-//template <class CONTAINER, class DATA>
void optimizeMesh(GRegion *gr, const qmTetrahedron::Measures &qm)
{
// well, this should not be true !!!
@@ -902,16 +901,9 @@ void optimizeMesh(GRegion *gr, const qmTetrahedron::Measures &qm)
// printf("coucou\n");
for (CONTAINER::iterator it = allTets.begin(); it != allTets.end(); ++it){
if (!(*it)->isDeleted()){
- for (int i=0;i<4;i++)
- for (int j=i+1;j<4;j++)
- if ((*it)->tet()->getVertex(i) == (*it)->tet()->getVertex(j))
- printf("argh\n");
-
double qq = (*it)->getQuality();
if (qq < qMin){
- // printf("cacze\n");
for (int i = 0; i < 6; i++){
- // printf("%d\n",i);
if (edgeSwap(newTets, *it, i, qm)) {
nbESwap++;
break;
diff --git a/Mesh/meshGRegionLocalMeshMod.cpp b/Mesh/meshGRegionLocalMeshMod.cpp
index 11e16d634957052530c1f400b2844a31ef8d1d9b..83488f67d1f3f7a66b8dfc9023010366ef1f2fd0 100644
--- a/Mesh/meshGRegionLocalMeshMod.cpp
+++ b/Mesh/meshGRegionLocalMeshMod.cpp
@@ -105,7 +105,10 @@ bool buildEdgeCavity(MTet4 *t, int iLocalEdge, MVertex **v1, MVertex **v2,
return false;
}
// FIXME when hybrid mesh, this loops for ever
- if (cavity.size() > 1000) return false;
+ if (cavity.size() > 1000) {
+ printf("cavity size gets laaaarge\n");
+ return false;
+ }
// printf("%d %d\n",ITER++, cavity.size());
}
computeNeighboringTetsOfACavity (cavity,outside);
diff --git a/Mesh/meshGRegionRelocateVertex.cpp b/Mesh/meshGRegionRelocateVertex.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1c004f84d768690e2268c2b481991ee2177b8d63
--- /dev/null
+++ b/Mesh/meshGRegionRelocateVertex.cpp
@@ -0,0 +1,119 @@
+#include "GModel.h"
+#include "GRegion.h"
+#include "MLine.h"
+#include "MTriangle.h"
+#include "MQuadrangle.h"
+#include "MTetrahedron.h"
+#include "MPyramid.h"
+#include "MPrism.h"
+#include "MHexahedron.h"
+#include "Context.h"
+#include "meshGFaceOptimize.h"
+static double objective_function (double xi, MVertex *ver,
+ double xTarget, double yTarget, double zTarget,
+ const std::vector<MElement*> <){
+ double x = ver->x();
+ double y = ver->y();
+ double z = ver->z();
+ ver->x() = (1.-xi) * ver->x() + xi * xTarget;
+ ver->y() = (1.-xi) * ver->y() + xi * yTarget;
+ ver->z() = (1.-xi) * ver->z() + xi * zTarget;
+ double minQual = 1.0;
+ for(unsigned int i = 0; i < lt.size(); i++){
+ if (lt[i]->getNumVertices () == 4)
+ minQual = std::min((lt[i]->minSICNShapeMeasure()), minQual);
+ else
+ minQual = std::min(fabs(lt[i]->minSICNShapeMeasure()), minQual);
+ }
+ ver->x() = x;
+ ver->y() = y;
+ ver->z() = z;
+ return minQual;
+}
+
+
+#define sqrt5 2.236067977499789696
+
+static int Stopping_Rule(double x0, double x1)
+{
+ double tolerance = 1.e-2;
+ return ( fabs( x1 - x0 ) < tolerance ) ? 1 : 0;
+}
+
+double Maximize_Quality_Golden_Section( MVertex *ver,
+ double xTarget, double yTarget, double zTarget,
+ const std::vector<MElement*> < )
+{
+
+ static const double lambda = 0.5 * (sqrt5 - 1.0);
+ static const double mu = 0.5 * (3.0 - sqrt5); // = 1 - lambda
+ double a = 0.0;
+ double b = 2.0;
+
+ double x1 = b - lambda * (b - a);
+ double x2 = a + lambda * (b - a);
+ double fx1 = objective_function (x1, ver, xTarget, yTarget, zTarget , lt );
+ double fx2 = objective_function (x2, ver, xTarget, yTarget, zTarget , lt );
+
+ while ( ! Stopping_Rule( a, b ) ) {
+ // printf("GOLDEN : %g %g (%12.5E,%12.5E)\n",a,b,fa,fb);
+ if (fx1 < fx2) {
+ a = x1;
+ if ( Stopping_Rule( a, b ) ) break;
+ x1 = x2;
+ fx1 = fx2;
+ x2 = b - mu * (b - a);
+ fx2 = objective_function (x2, ver, xTarget, yTarget, zTarget , lt );
+ } else {
+ b = x2;
+ if ( Stopping_Rule( a, b ) ) break;
+ x2 = x1;
+ fx2 = fx1;
+ x1 = a + mu * (b - a);
+ fx1 = objective_function (x1, ver, xTarget, yTarget, zTarget , lt );
+ }
+ }
+ // printf("finally : %g %g (%12.5E,%12.5E)\n",a,b,fa,fb);
+ return a;
+}
+
+static void _relocateVertex(MVertex *ver,
+ const std::vector<MElement*> <)
+{
+ if(ver->onWhat()->dim() != 3) return;
+ double x = 0, y=0, z=0;
+ int N = 0;
+ for(unsigned int i = 0; i < lt.size(); i++){
+ double XCG=0,YCG=0,ZCG=0;
+ for (int j=0;j<lt[i]->getNumVertices();j++){
+ XCG += lt[i]->getVertex(j)->x();
+ YCG += lt[i]->getVertex(j)->y();
+ ZCG += lt[i]->getVertex(j)->z();
+ }
+ x += XCG;
+ y += YCG;
+ z += ZCG;
+ N += lt[i]->getNumVertices();
+ }
+
+ double xi = Maximize_Quality_Golden_Section( ver, x/N, y/N, z/N, lt );
+ ver->x() = (1.-xi) * ver->x() + xi * x/N;
+ ver->y() = (1.-xi) * ver->y() + xi * y/N;
+ ver->z() = (1.-xi) * ver->z() + xi * z/N;
+}
+
+void RelocateVertices (std::vector<GRegion*> ®ions) {
+ for(unsigned int k = 0; k < regions.size(); k++){
+ v2t_cont adj;
+ buildVertexToElement(regions[k]->tetrahedra, adj);
+ buildVertexToElement(regions[k]->pyramids, adj);
+ buildVertexToElement(regions[k]->prisms, adj);
+ buildVertexToElement(regions[k]->hexahedra, adj);
+ v2t_cont::iterator it = adj.begin();
+ while (it != adj.end()){
+ _relocateVertex( it->first, it->second);
+ ++it;
+ }
+ }
+}
+
diff --git a/Mesh/meshGRegionRelocateVertex.h b/Mesh/meshGRegionRelocateVertex.h
new file mode 100644
index 0000000000000000000000000000000000000000..678403aade1885bc07aa04f336f5d0c5e9e2a1dd
--- /dev/null
+++ b/Mesh/meshGRegionRelocateVertex.h
@@ -0,0 +1,6 @@
+#ifndef _MESHGREGIONRELOCATEVERTEX_
+#define _MESHGREGIONRELOCATEVERTEX_
+#include <vector>
+class GRegion;
+void RelocateVertices (std::vector<GRegion*> ®ions);
+#endif
diff --git a/benchmarks/3d/hex.geo b/benchmarks/3d/hex.geo
index 8b533e4d68cdacb5b0f1db00ed68cb1ffe330a47..f6ff4b542e49c20ba043490541e39a93c300c467 100644
--- a/benchmarks/3d/hex.geo
+++ b/benchmarks/3d/hex.geo
@@ -1,4 +1,4 @@
-Mesh.CharacteristicLengthExtendFromBoundary = 0;
+//Mesh.CharacteristicLengthExtendFromBoundary = 0;
lc = 0.3;
// example of a purely hexahedral mesh using only transfinite
@@ -40,17 +40,17 @@ Line Loop(24) = {6,-12,3,10};
Ruled Surface(25) = {24};
Surface Loop(1) = {17,-25,-23,-21,19,15};
Volume(1) = {1};
-Transfinite Line{1:4,6:9} = 30;
-Transfinite Line{10:13} = 30;
+//Transfinite Line{1:4,6:9} = 20;
+//Transfinite Line{10:13} = 20;
-Transfinite Surface {15} = {1,2,3,4};
+//Transfinite Surface {15} = {1,2,3,4};
-Transfinite Surface {17} = {5,6,7,8};
-Transfinite Surface {19} = {1,5,8,4};
-Transfinite Surface {21} = {8,7,3,4};
-Transfinite Surface {23} = {6,7,3,2};
-Transfinite Surface {25} = {5,6,2,1};
+//Transfinite Surface {17} = {5,6,7,8};
+//Transfinite Surface {19} = {1,5,8,4};
+//Transfinite Surface {21} = {8,7,3,4};
+//Transfinite Surface {23} = {6,7,3,2};
+//Transfinite Surface {25} = {5,6,2,1};
//Transfinite Volume{1} = {1,2,3,4,5,6,7,8};