// Gmsh - Copyright (C) 1997-2017 C. Geuzaine, J.-F. Remacle // // See the LICENSE.txt file for license information. Please report all // bugs and problems to the public mailing list <gmsh@onelab.info>. #include "boundaryLayersData.h" #include "GmshConfig.h" #include "GModel.h" #include "MLine.h" #include "MTriangle.h" #include "MEdge.h" #include "OS.h" #if !defined(HAVE_MESH) || !defined(HAVE_ANN) BoundaryLayerField* getBLField(GModel *gm){ return 0; } bool buildAdditionalPoints2D (GFace *gf ) { return false; } bool buildAdditionalPoints3D (GRegion *gr) { return false; } //void buildMeshMetric(GFace *gf, double *uv, SMetric3 &m, double metric[3]) {} /*faceColumn BoundaryLayerColumns::getColumns(GFace *gf, MVertex *v1, MVertex *v2, MVertex *v3, int side) { return faceColumn(BoundaryLayerData(),BoundaryLayerData(),BoundaryLayerData()); } */ edgeColumn BoundaryLayerColumns::getColumns(MVertex *v1, MVertex *v2 , int side) { return edgeColumn(BoundaryLayerData(),BoundaryLayerData()); } #else #include "Field.h" const int FANSIZE__ = 4; /* ^ ni | | +-----------------+ bi / bj / /\ / \ nj / Z + */ SVector3 interiorNormal(SPoint2 p1, SPoint2 p2, SPoint2 p3) { SVector3 ez (0,0,1); SVector3 d (p1.x()-p2.x(),p1.y()-p2.y(),0); SVector3 h (p3.x()-0.5*(p2.x()+p1.x()),p3.y()-0.5*(p2.y()+p1.y()),0); SVector3 n = crossprod(d,ez); n.normalize(); if (dot(n,h) > 0)return n; return n*(-1.); } edgeColumn BoundaryLayerColumns::getColumns(MVertex *v1, MVertex *v2 , int side) { Equal_Edge aaa; MEdge e(v1, v2); std::map<MVertex*,BoundaryLayerFan>::const_iterator it1 = _fans.find(v1); std::map<MVertex*,BoundaryLayerFan>::const_iterator it2 = _fans.find(v2); int N1 = getNbColumns(v1) ; int N2 = getNbColumns(v2) ; int nbSides = _normals.count(e); // if (nbSides != 1)printf("I'm here %d sides\n",nbSides); // Standard case, only two extruded columns from the two vertices if (N1 == 1 && N2 == 1) return edgeColumn(getColumn(v1,0),getColumn(v2,0)); // one fan on if (nbSides == 1){ if (it1 != _fans.end() && it2 == _fans.end() ){ if (aaa(it1->second._e1,e)) return edgeColumn(getColumn (v1,0),getColumn(v2,0)); else return edgeColumn(getColumn (v1,N1-1),getColumn(v2,0)); } if (it2 != _fans.end() && it1 == _fans.end() ){ if (aaa(it2->second._e1,e)) return edgeColumn(getColumn (v1,0),getColumn(v2,0)); else return edgeColumn(getColumn (v1,0),getColumn(v2,N2-1)); } if (it2 != _fans.end() && it1 != _fans.end() ){ int c1, c2; if (aaa(it1->second._e1,e)) c1 = 0; else c1 = N1-1; if (aaa(it2->second._e1,e)) c2 = 0; else c2 = N2-1; return edgeColumn(getColumn (v1,c1),getColumn(v2,c2)); } // fan on the right if (N1 == 1 || N2 == 2){ const BoundaryLayerData & c10 = getColumn(v1,0); const BoundaryLayerData & c20 = getColumn(v2,0); const BoundaryLayerData & c21 = getColumn(v2,1); if (dot(c10._n,c20._n) > dot(c10._n,c21._n) ) return edgeColumn(c10,c20); else return edgeColumn(c10,c21); } // fan on the left if (N1 == 2 || N2 == 1){ const BoundaryLayerData & c10 = getColumn(v1,0); const BoundaryLayerData & c11 = getColumn(v1,1); const BoundaryLayerData & c20 = getColumn(v2,0); if (dot(c10._n,c20._n) > dot(c11._n,c20._n) ) return edgeColumn(c10,c20); else return edgeColumn(c11,c20); } // Msg::Error ("Impossible Boundary Layer Configuration : " // "one side and no fans %d %d", N1, N2); // FIXME WRONG return N1 ? edgeColumn (getColumn (v1,0),getColumn(v1,0)) : edgeColumn (getColumn (v2,0),getColumn(v2,0)); } else if (nbSides == 2){ int i1=0,i2=1,j1=0,j2=1; if (it1 != _fans.end()){ i1 = aaa(it1->second._e1,e) ? 0 : N1-1; i2 = !aaa(it1->second._e1,e) ? 0 : N1-1; } if (it2 != _fans.end()){ j1 = aaa(it2->second._e1,e) ? 0 : N2-1; j2 = !aaa(it2->second._e1,e) ? 0 : N2-1; } const BoundaryLayerData & c10 = getColumn(v1,i1); const BoundaryLayerData & c11 = getColumn(v1,i2); const BoundaryLayerData & c20 = getColumn(v2,j1); const BoundaryLayerData & c21 = getColumn(v2,j2); if (side == 0){ if (dot(c10._n,c20._n) > dot(c10._n,c21._n) ) return edgeColumn(c10,c20); else return edgeColumn(c10,c21); } if (side == 1){ if (dot(c11._n,c20._n) > dot(c11._n,c21._n) ) return edgeColumn(c11,c20); else return edgeColumn(c11,c21); } } Msg::Error("Not yet Done in BoundaryLayerData nbSides = %d, ",nbSides ); static BoundaryLayerData error; static edgeColumn error2(error, error); return error2; } /* static bool pointInFace (GFace *gf, double u, double v) { return true; } */ static void treat2Connections(GFace *gf, MVertex *_myVert, MEdge &e1, MEdge &e2, BoundaryLayerColumns *_columns, std::vector<SVector3> &_dirs, bool fan) { std::vector<SVector3> N1,N2; for (std::multimap<MEdge,SVector3,Less_Edge>::iterator itm = _columns->_normals.lower_bound(e1); itm != _columns->_normals.upper_bound(e1); ++itm) N1.push_back(itm->second); for (std::multimap<MEdge,SVector3,Less_Edge>::iterator itm = _columns->_normals.lower_bound(e2); itm != _columns->_normals.upper_bound(e2); ++itm) N2.push_back(itm->second); if (N1.size() == N2.size()){ for (unsigned int SIDE = 0; SIDE < N1.size() ; SIDE++){ if (!fan){ SVector3 x = N1[SIDE]*1.01+N2[SIDE]; x.normalize(); _dirs.push_back(x); } else if (fan){ // printf("fan \n"); int fanSize = FANSIZE__; // if the angle is greater than PI, than reverse the sense double alpha1 = atan2(N1[SIDE].y(),N1[SIDE].x()); double alpha2 = atan2(N2[SIDE].y(),N2[SIDE].x()); double AMAX = std::max(alpha1,alpha2); double AMIN = std::min(alpha1,alpha2); MEdge ee[2]; if (alpha1 > alpha2){ ee[0] = e2;ee[1] = e1; } else { ee[0] = e1;ee[1] = e2; } if ( AMAX - AMIN >= M_PI){ double temp = AMAX; AMAX = AMIN + 2*M_PI; AMIN = temp; MEdge eee0 = ee[0]; ee[0] = ee[1];ee[1] = eee0; } _columns->addFan (_myVert,ee[0],ee[1],true); for (int i=-1; i<=fanSize; i++){ double t = (double)(i+1)/ (fanSize+1); double alpha = t * AMAX + (1.-t)* AMIN; SVector3 x (cos(alpha),sin(alpha),0); x.normalize(); _dirs.push_back(x); } } /* else { _dirs.push_back(N1[SIDE]); _dirs.push_back(N2[SIDE]); } */ } } } static void treat3Connections(GFace *gf, MVertex *_myVert, MEdge &e1, MEdge &e2, MEdge &e3, BoundaryLayerColumns *_columns, std::vector<SVector3> &_dirs) { std::vector<SVector3> N1,N2,N3; for (std::multimap<MEdge,SVector3,Less_Edge>::iterator itm = _columns->_normals.lower_bound(e1); itm != _columns->_normals.upper_bound(e1); ++itm) N1.push_back(itm->second); for (std::multimap<MEdge,SVector3,Less_Edge>::iterator itm = _columns->_normals.lower_bound(e2); itm != _columns->_normals.upper_bound(e2); ++itm) N2.push_back(itm->second); for (std::multimap<MEdge,SVector3,Less_Edge>::iterator itm = _columns->_normals.lower_bound(e3); itm != _columns->_normals.upper_bound(e3); ++itm) N3.push_back(itm->second); SVector3 x1,x2; if (N1.size() == 2){ } else if (N2.size() == 2){ std::vector<SVector3> temp = N1; N1.clear(); N1 = N2; N2.clear(); N2 = temp; } else if (N3.size() == 2){ std::vector<SVector3> temp = N1; N1.clear(); N1 = N3; N3.clear(); N3 = temp; } else { Msg::Fatal("IMPOSSIBLE BL CONFIGURATION"); } if (dot(N1[0],N2[0]) > dot(N1[0],N3[0])){ x1 = N1[0]*1.01+N2[0]; x2 = N1[1]*1.01+N3[0]; } else { x1 = N1[1]*1.01+N2[0]; x2 = N1[0]*1.01+N3[0]; } x1.normalize(); _dirs.push_back(x1); x2.normalize(); _dirs.push_back(x2); } BoundaryLayerField *getBLField(GModel *gm) { FieldManager *fields = gm->getFields(); if(fields->getBoundaryLayerField() <= 0){ return 0; } Field *bl_field = fields->get(fields->getBoundaryLayerField()); return dynamic_cast<BoundaryLayerField*> (bl_field); } static bool isEdgeOfFaceBL(GFace *gf, GEdge *ge, BoundaryLayerField *blf) { if (blf->isEdgeBL (ge->tag()))return true; /* std::list<GFace*> faces = ge->faces(); for (std::list<GFace*>::iterator it = faces.begin(); it != faces.end() ; ++it){ if ((*it) == gf)return false; } for (std::list<GFace*>::iterator it = faces.begin(); it != faces.end() ; ++it){ if (blf->isFaceBL ((*it)->tag()))return true; } */ return false; } static void getEdgesData(GFace *gf, BoundaryLayerField *blf, BoundaryLayerColumns *_columns, std::set<MVertex*> &_vertices, std::set<MEdge,Less_Edge> &allEdges, std::multimap<MVertex*,MVertex*> &tangents) { // get all model edges std::list<GEdge*> edges = gf->edges(); std::list<GEdge*> embedded_edges = gf->embeddedEdges(); edges.insert(edges.begin(), embedded_edges.begin(),embedded_edges.end()); // iterate on model edges std::list<GEdge*>::iterator ite = edges.begin(); while(ite != edges.end()){ // check if this edge generates a boundary layer if (isEdgeOfFaceBL (gf,*ite,blf)){ for(unsigned int i = 0; i< (*ite)->lines.size(); i++){ MVertex *v1 = (*ite)->lines[i]->getVertex(0); MVertex *v2 = (*ite)->lines[i]->getVertex(1); allEdges.insert(MEdge(v1,v2)); _columns->_non_manifold_edges.insert(std::make_pair(v1,v2)); _columns->_non_manifold_edges.insert(std::make_pair(v2,v1)); _vertices.insert(v1); _vertices.insert(v2); } } else { MVertex *v1 = (*ite)->lines[0]->getVertex(0); MVertex *v2 = (*ite)->lines[0]->getVertex(1); MVertex *v3 = (*ite)->lines[(*ite)->lines.size() - 1]->getVertex(1); MVertex *v4 = (*ite)->lines[(*ite)->lines.size() - 1]->getVertex(0); tangents.insert (std::make_pair (v1,v2)); tangents.insert (std::make_pair (v3,v4)); } ++ite; } } static void getNormals(GFace *gf, BoundaryLayerField *blf, BoundaryLayerColumns *_columns, std::set<MEdge,Less_Edge> &allEdges) { // assume that the initial mesh has been created i.e. that there exist // triangles inside the domain. Triangles are used to define // exterior normals for (unsigned int i = 0; i < gf->triangles.size(); i++){ SPoint2 p0,p1,p2; MVertex *v0 = gf->triangles[i]->getVertex(0); MVertex *v1 = gf->triangles[i]->getVertex(1); MVertex *v2 = gf->triangles[i]->getVertex(2); reparamMeshEdgeOnFace(v0, v1, gf, p0, p1); reparamMeshEdgeOnFace(v0, v2, gf, p0, p2); MEdge me01(v0,v1); if (allEdges.find(me01) != allEdges.end()){ SVector3 v01 = interiorNormal (p0,p1,p2); _columns->_normals.insert(std::make_pair(me01,v01)); } MEdge me02(v0,v2); if (allEdges.find(me02) != allEdges.end()){ SVector3 v02 = interiorNormal (p0,p2,p1); _columns->_normals.insert(std::make_pair(me02,v02)); } MEdge me21(v2,v1); if (allEdges.find(me21) != allEdges.end()){ SVector3 v21 = interiorNormal (p2,p1,p0); _columns->_normals.insert(std::make_pair(me21,v21)); } } } static void addColumnAtTheEndOfTheBL(GEdge *ge, GVertex *gv, BoundaryLayerColumns* _columns, BoundaryLayerField *blf) { if (!blf->isEdgeBL(ge->tag())){ GVertex *g0 = ge->getBeginVertex(); GVertex *g1 = ge->getEndVertex(); MVertex * v0 = g0->mesh_vertices[0]; MVertex * v1 = g1->mesh_vertices[0]; std::vector<MVertex*> invert; for(unsigned int i = 0; i < ge->mesh_vertices.size() ; i++) invert.push_back(ge->mesh_vertices[ge->mesh_vertices.size() - i - 1]); SVector3 t (v1->x()-v0->x(), v1->y()-v0->y(),v1->z()-v0->z()); t.normalize(); if (g0 == gv) _columns->addColumn(t, v0, ge->mesh_vertices); else if (g1 == gv) _columns->addColumn(t*-1.0, v1,invert); } } void getLocalInfoAtNode (MVertex *v, BoundaryLayerField *blf, double &hwall) { hwall = blf->hwall_n; if (v->onWhat()->dim() == 0){ hwall= blf->hwall (v->onWhat()->tag()); } else if (v->onWhat()->dim() == 1){ GEdge *ge = (GEdge*)v->onWhat(); Range<double> bounds = ge->parBounds(0); double t_begin = bounds.low(); double t_end = bounds.high(); double t; v->getParameter(0,t); double hwall_beg = blf->hwall (ge->getBeginVertex()->tag()); double hwall_end = blf->hwall (ge->getEndVertex()->tag()); hwall = hwall_beg + (t-t_begin)/(t_end-t_begin) * (hwall_end - hwall_beg); } } bool buildAdditionalPoints2D(GFace *gf) { BoundaryLayerColumns *_columns = gf->getColumns(); _columns->clearData(); // GET THE FIELD THAT DEFINES THE DISTANCE FUNCTION BoundaryLayerField *blf = getBLField (gf->model()); if (!blf)return false; blf->setupFor2d(gf->tag()); std::set<MVertex*> _vertices; std::set<MEdge,Less_Edge> allEdges; std::multimap<MVertex*,MVertex*> tangents; getEdgesData ( gf, blf, _columns, _vertices , allEdges , tangents ); if (!_vertices.size())return false; getNormals ( gf, blf, _columns, allEdges); // for all boundry points for (std::set<MVertex*>::iterator it = _vertices.begin(); it != _vertices.end() ; ++it){ bool endOfTheBL = false; SVector3 dirEndOfBL; std::vector<MVertex*> columnEndOfBL; std::vector<MVertex*> _connections; std::vector<SVector3> _dirs; // get all vertices that are connected to that // vertex among all boundary layer vertices ! bool fan = (*it)->onWhat()->dim() == 0 && blf->isFanNode((*it)->onWhat()->tag()); for (std::multimap<MVertex*,MVertex*>::iterator itm = _columns->_non_manifold_edges.lower_bound(*it); itm != _columns->_non_manifold_edges.upper_bound(*it); ++itm) _connections.push_back (itm->second); // A trailing edge topology : 3 edges incident to a vertex if (_connections.size() == 3){ MEdge e1 (*it,_connections[0]); MEdge e2 (*it,_connections[1]); MEdge e3 (*it,_connections[2]); treat3Connections (gf, *it,e1,e2,e3, _columns, _dirs); } // STANDARD CASE, one vertex connected to two neighboring vertices else if (_connections.size() == 2){ MEdge e1 (*it,_connections[0]); MEdge e2 (*it,_connections[1]); treat2Connections (gf, *it,e1,e2, _columns, _dirs, fan); } else if (_connections.size() == 1){ MEdge e1 (*it,_connections[0]); std::vector<SVector3> N1; for (std::multimap<MEdge,SVector3,Less_Edge>::iterator itm = _columns->_normals.lower_bound(e1); itm != _columns->_normals.upper_bound(e1); ++itm) N1.push_back(itm->second); // one point has only one side and one normal : it has to be at the end of the BL // then, we have the tangent to the connecting edge // *it _connections[0] // --------- + ----------- // NO BL BL if (N1.size() == 1){ std::vector<MVertex*> Ts; for (std::multimap<MVertex*,MVertex*>::iterator itm = tangents.lower_bound(*it); itm != tangents.upper_bound(*it); ++itm) Ts.push_back(itm->second); // end of the BL --> let's add a column that correspond to the // model edge that lies after the end of teh BL if (Ts.size() == 1){ // printf("HERE WE ARE IN FACE %d %d\n",gf->tag(),Ts.size()); // printf("Classif dim %d %d\n",(*it)->onWhat()->dim(),Ts[0]->onWhat()->dim()); GEdge *ge = dynamic_cast<GEdge*>(Ts[0]->onWhat()); GVertex *gv = dynamic_cast<GVertex*>((*it)->onWhat()); if (ge && gv){ addColumnAtTheEndOfTheBL (ge,gv,_columns,blf); } } else { Msg::Error("Impossible BL Configuration -- One Edge -- Tscp.size() = %d",Ts.size()); } } else if (N1.size() == 2){ // printf("%g %g --> %g %g \n",e1.getVertex(0)->x(),e1.getVertex(0)->y(), // e1.getVertex(1)->x(),e1.getVertex(1)->y()); // printf("N1.size = %d %g %g %g %g\n",N1.size(),N1[0].x(),N1[0].y(),N1[1].x(),N1[1].y()); SPoint2 p0,p1; reparamMeshEdgeOnFace(*it,_connections[0], gf, p0, p1); int fanSize = FANSIZE__; double alpha1 = atan2(N1[0].y(),N1[0].x()); double alpha2 = atan2(N1[1].y(),N1[1].x()); double alpha3 = atan2(p1.y()-p0.y(),p1.x()-p0.x()); double AMAX = std::max(alpha1,alpha2); double AMIN = std::min(alpha1,alpha2); if (alpha3 > AMAX){ AMIN += M_PI; AMAX += M_PI; } if ( AMAX - AMIN >= M_PI){ double temp = AMAX; AMAX = AMIN + 2*M_PI; AMIN = temp; } _columns->addFan (*it,e1,e1,true); // printf("%g %g --> %g %g\n",N1[0].x(),N1[0].y(),N1[1].x(),N1[1].y()); for (int i=-1; i<=fanSize; i++){ double t = (double)(i+1)/ (fanSize+1); double alpha = t * AMAX + (1.-t)* AMIN; SVector3 x (cos(alpha),sin(alpha),0); // printf("%d %g %g %g\n",i,x.x(),x.y(),alpha); x.normalize(); _dirs.push_back(x); } } } // if (_dirs.size() > 1)printf("%d directions\n",_dirs.size()); // now create the BL points for (unsigned int DIR=0;DIR<_dirs.size();DIR++){ SVector3 n = _dirs[DIR]; // < ------------------------------- > // // N = X(p0+ e n) - X(p0) // // = e * (dX/du n_u + dX/dv n_v) // // < ------------------------------- > // /* if (endOfTheBL){ printf("%g %g %d %d %g\n", (*it)->x(), (*it)->y(), DIR, (int)_dirs.size(), dot(n, dirEndOfBL)); } */ if (endOfTheBL && dot(n,dirEndOfBL) > .99){ // printf( "coucou c'est moi\n"); } else { MVertex *first = *it; double hwall; getLocalInfoAtNode (first, blf, hwall); std::vector<MVertex*> _column; SPoint2 par = gf->parFromPoint(SPoint3(first->x(),first->y(),first->z())); double L = hwall; while(1){ // printf("L = %g\n",L); if (L > blf->thickness) break; SPoint2 pnew (par.x() + L * n.x(), par.y() + L * n.y()); GPoint pp = gf->point(pnew); MFaceVertex *_current = new MFaceVertex (pp.x(),pp.y(),pp.z(),gf,pnew.x(),pnew.y()); _current->bl_data = new MVertexBoundaryLayerData; _column.push_back(_current); int ith = _column.size() ; L+= hwall * pow (blf->ratio, ith); } _columns->addColumn(n,*it, _column /*,_metrics*/); } } } // DEBUG STUFF char name[256]; sprintf(name, "points_face_%d.pos", gf->tag()); FILE *f = Fopen (name,"w"); if(f){ fprintf(f,"View \"\" {\n"); for (std::set<MVertex*>::iterator it = _vertices.begin(); it != _vertices.end() ; ++it){ MVertex *v = *it; for (int i=0;i<_columns->getNbColumns(v);i++){ const BoundaryLayerData &data = _columns->getColumn(v,i); for (unsigned int j = 0; j < data._column.size(); j++){ MVertex *blv = data._column[j]; fprintf(f,"SP(%g,%g,%g){%d};\n",blv->x(),blv->y(),blv->z(),v->getNum()); } } } fprintf(f,"};\n"); fclose (f); } // END OF DEBUG STUFF return 1; } #endif