1 /*
2 XLiFE++ is an extended library of finite elements written in C++
3 Copyright (C) 2014 Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin
4
5 This program is free software: you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13 You should have received a copy of the GNU General Public License
14 along with this program. If not, see <http://www.gnu.org/licenses/>.
15 */
16
17 /*!
18 \file GeomFigureMesh.cpp
19 \author Y. Lafranche
20 \since 22 May 2008
21 \date 08 Mar 2014
22
23 \brief Implementation of xlifepp::subdivision::GeomFigureMesh class members and related functions
24 */
25
26 #include "GeomFigureMesh.hpp"
27 #include "Hexahedron.hpp"
28 #include "Tetrahedron.hpp"
29 #include "Quadrangle.hpp"
30 #include "Triangle.hpp"
31 #include "TeXPolygon.hpp"
32
33 #include <algorithm>
34 #include<vector>
35 using namespace std;
36
37 namespace xlifepp {
38 namespace subdivision {
39
40 //-------------------------------------------------------------------------------
41 // Public member functions
42 //-------------------------------------------------------------------------------
43
44 /*!
45 Returns the list of all the vertices of the element num, with
46 minElementNum_ <= num <= numberOfElements+minElementNum_-1.
47 Each vertex is given by its number, which is >= minVertexNum_.
48 */
49 template<class T_>
element(const number_t num) const50 vector<number_t> GeomFigureMesh<T_>::element(const number_t num) const {
51 vector<number_t> V = listT_.at(num-minElementNum_).rankOfVertices();
52 rankToNum(V);
53 return V;
54 }
55
56 /*!
57 returns the localization code of the element num, with num such that
58 minElementNum_ <= num <= numberOfElements+minElementNum_-1
59 The code is computed from those of its main vertices. The function
60 returns the localization code of the domain the element belongs to.
61 */
62 template<class T_>
locCode(const number_t num) const63 refnum_t GeomFigureMesh<T_>::locCode(const number_t num) const {
64 const T_& elem(listT_.at(num-minElementNum_));
65 vector<number_t> V = elem.rankOfVertices();
66 refnum_t localcod = listV_[V[0]].locCode();
67 for (number_t j=1; j<elem.numberOfO1Vertices(); j++) {
68 localcod &= listV_[V[j]].locCode();
69 }
70 return localcod;
71 }
72
73 /*!
74 returns a normal vector to the face number i >= 1 of element num,
75 with num such that minElementNum_ <= num <= numberOfElements+minElementNum_-1
76 The vector is not normalized and is oriented towards the exterior of the element.
77 */
78 template<class T_>
faceExtNormalVec(const number_t num,const number_t i) const79 vector<real_t> GeomFigureMesh<T_>::faceExtNormalVec(const number_t num, const number_t i) const {
80 return listT_.at(num-minElementNum_).extNormVec(i, listV_);
81 }
82
83 /*!
84 returns the orientation of the face number i >= 1 of element num, with num such that
85 minElementNum_ <= num <= numberOfElements+minElementNum_-1. The function:
86 - returns 1 if the first three vertices (i,j,k) of the face define a vector ij x ik
87 oriented towards the exterior of the element,
88 - returns -1 otherwise.
89 */
90 template<class T_>
faceOrientation(const number_t num,const number_t i) const91 int GeomFigureMesh<T_>::faceOrientation(const number_t num, const number_t i) const {
92 return listT_.at(num-minElementNum_).faceOrientation(i);
93 }
94
95 /*!
96 Returns the list of elements of the mesh belonging to an area of kind TA
97 or belonging to any area of kind TA if num=0. The elements are given by their
98 number in the whole set of elements of the mesh.
99 */
100 template<class T_>
elementsIn(const topologicalArea TA,const number_t num) const101 vector<number_t> GeomFigureMesh<T_>::elementsIn(const topologicalArea TA, const number_t num) const {
102 vector<number_t> V;
103 refnum_t sig = lCodeOf(TA,num);
104
105 typename vector<T_>::const_iterator itT;
106 for (itT=listT_.begin(); itT != listT_.end(); itT++) {
107 refnum_t localcod = listV_[itT->rankOfVertex(1)].locCode();
108 for (number_t j=2; j<=nb_main_vertices_by_element_; j++) {
109 localcod &= listV_[itT->rankOfVertex(j)].locCode();
110 }
111 // If the element belongs to the area, store its number.
112 if (localcod & sig) {V.push_back(itT->number());}
113 }
114 return V;
115 }
116
117 /*!
118 Returns the number of edges belonging to an area num of kind TA
119 or belonging to any area of this kind if num=0.
120 */
121 template<class T_>
numberOfEdgesIn(const topologicalArea TA,const number_t num) const122 number_t GeomFigureMesh<T_>::numberOfEdgesIn(const topologicalArea TA, const number_t num) const {
123 return rk_edgesIn(TA,num).size();
124 }
125 /*!
126 Returns the list of edges of the elements of the mesh belonging to an area num
127 of kind TA or belonging to any area of this kind if num=0.
128 Each edge is defined by a pair of the form (nelem[i], nedge[i]) where
129 . nelem[i] is the number of the element (>=minElementNum_) the edge belongs to.
130 . nedge[i] is the local number of the edge in the element (from 1 to nb_edges_by_element_)
131 These numbers are stored in two separate vectors nelem and nedge, with the same
132 length, grouped as a pair returned by the function.
133 -> The core of the algorithm is the same in rk_edgesIn.
134 Nota: If TA is an interface, the result consists in fact in two lists, since each
135 segment of the interface belongs to 2 elements. The two lists are intermixed,
136 except when there are only two subdomains, in which case, due to the construction
137 algorithm, the first half of the elements belong to one subdomain and the second
138 half of the elements belong to the other subdomain.
139 If needed, the two lists should be separated using geometric computations.
140 */
141 template<class T_>
142 pair< vector<number_t>,vector<number_t> >
edgeElementsIn(const topologicalArea TA,const number_t num) const143 GeomFigureMesh<T_>::edgeElementsIn(const topologicalArea TA, const number_t num) const {
144 vector<number_t> numEdge;
145 vector<number_t> numElem;
146 refnum_t sig = lCodeOf(TA,num);
147
148 typename vector<T_>::const_iterator itT;
149 for (itT=listT_.begin(); itT != listT_.end(); itT++) {
150 for (number_t i=1; i<=nb_edges_by_element_; i++) {
151 pair_nn prV = itT->rkOfO1VeOnEdge(i);
152 refnum_t localcod = listV_[prV.first].locCode() & listV_[prV.second].locCode();
153 // If the edge belongs to the area, store it.
154 if (localcod & sig) { numEdge.push_back(i); numElem.push_back(itT->number()); }
155 }
156 }
157 return make_pair(numElem, numEdge);
158 }
159 /*!
160 Returns the list of edges of a triangle mesh belonging to an area num of kind TA
161 or belonging to any area of this kind if num=0.
162 Edges are defined by the number of their vertices (starting from minVertexNum_).
163 */
164 template<class T_>
edgesIn(const topologicalArea TA,const number_t num) const165 vector< pair_nn > GeomFigureMesh<T_>::edgesIn(const topologicalArea TA, const number_t num) const {
166 vector< pair_nn > EonB = rk_edgesIn(TA,num);
167 vector< pair_nn >::iterator itEonB;
168
169 for (itEonB=EonB.begin(); itEonB != EonB.end(); itEonB++) { rankToNum(*itEonB); }
170 return EonB;
171 }
172
173 // functions for 3D geometries
174 /*!
175 Returns the number of faces belonging to an area num of kind TA
176 or belonging to any area of this kind if num=0.
177 */
178 template<class T_>
numberOfFacesIn(const topologicalArea TA,const number_t num) const179 number_t GeomFigureMesh<T_>::numberOfFacesIn(const topologicalArea TA, const number_t num) const {
180 return rk_facesIn(TA,num).size();
181 }
182 /*!
183 Returns the number of the internal faces of the mesh
184 */
185 template<class T_>
numberOfFacesInside() const186 number_t GeomFigureMesh<T_>::numberOfFacesInside() const {
187 // lazy algorithm: it would be more efficient not to buid the list in memory,
188 // but just count the internal faces as they appear in the loop
189 // (see function facesInside).
190 return facesInside().size();
191 }
192
193 /*!
194 Returns the list of faces of the elements of the mesh belonging to an area num
195 of kind TA or belonging to any area of this kind if num=0.
196 Each face is defined by a pair of the form (nelem[i], nface[i]) where
197 . nelem[i] is the number of the element (>=minElementNum_) the face belongs to,
198 . nface[i] is the local number of the face in the element (from 1 to nb_faces_by_element_).
199 These numbers are stored in two separate vectors nelem and nface, with the same
200 length, grouped as a pair returned by the function.
201 -> The core of the algorithm is the same in rk_facesIn.
202 Nota: If TA is an interface, the result consists in fact in two lists, since each 2D
203 element of the interface belongs to two 3D elements. The two lists are intermixed,
204 except when there are only two subdomains, in which case, due to the construction
205 algorithm, the first half of the elements belong to one subdomain and the second
206 half of the elements belong to the other subdomain.
207 If needed, the two lists should be separated using geometric computations.
208 */
209 template<class T_>
210 pair< vector<number_t>,vector<number_t> >
faceElementsIn(const topologicalArea TA,const number_t num) const211 GeomFigureMesh<T_>::faceElementsIn(const topologicalArea TA, const number_t num) const {
212 vector<number_t> numFace;
213 vector<number_t> numElem;
214 refnum_t sig = lCodeOf(TA,num);
215
216 typename vector<T_>::const_iterator itT;
217 for (itT=listT_.begin(); itT != listT_.end(); itT++) {
218 for (number_t i=1; i<=nb_faces_by_element_; i++) {
219 vector<number_t> vrV = itT->rkOfO1VeOnFace(i);
220 vector<number_t>::const_iterator itvrV(vrV.begin());
221 refnum_t localcod = listV_[*itvrV].locCode();
222 while (++itvrV < vrV.end()) { localcod &= listV_[*itvrV].locCode(); }
223 // If the face belongs to the area, store it.
224 if (localcod & sig) { numFace.push_back(i); numElem.push_back(itT->number()); }
225 }
226 }
227 return make_pair(numElem, numFace);
228 }
229
230 /*!
231 Returns the list of faces of a mesh belonging to an area num of kind TA
232 or belonging to any area of this kind if num=0.
233 Faces are defined by the number of their vertices (starting from minVertexNum_).
234 */
235 template<class T_>
facesIn(const topologicalArea TA,const number_t num) const236 vector< vector<number_t> > GeomFigureMesh<T_>::facesIn(const topologicalArea TA, const number_t num) const {
237 vector< vector<number_t> > FonB = rk_facesIn(TA,num);
238 vector< vector<number_t> >::iterator itFonB;
239
240 for (itFonB=FonB.begin(); itFonB != FonB.end(); itFonB++) { rankToNum(*itFonB); }
241 return FonB;
242 }
243
244 /*!
245 Returns the list of faces of a mesh which are inside the boundaries.
246 Each face is defined by a pair of the form ( vertices, VnTnF ) where
247 . vertices is a set of integers, corresponding to the numbers of the vertices
248 defining the face,
249 . VnTnF is a vector of two pairs of the form (nelem, nface) corresponding to the
250 two elements sharing the face, where
251 . nelem is the number (>=1) of one of the two elements the face belongs to,
252 . nface is the local number of the face in this element (from 1 to nb_faces_by_element_).
253 */
254 template<class T_>
facesInside() const255 map_set_vec GeomFigureMesh<T_>::facesInside() const {
256 map_set_vec FI;
257 pair< map_set_vec::iterator, bool > retval;
258 refnum_t BMask = TG_.maskOf(boundaryArea);
259
260 typename vector<T_>::const_iterator itT;
261 for (itT=listT_.begin(); itT != listT_.end(); itT++) {
262 for (number_t i=1; i<=nb_faces_by_element_; i++) {
263 vector<number_t> vrV = itT->rkOfO1VeOnFace(i);
264 vector<number_t>::const_iterator itvrV(vrV.begin());
265 refnum_t localcod = BMask & listV_[*itvrV].locCode();
266 while (++itvrV < vrV.end()) { localcod &= listV_[*itvrV].locCode(); }
267 // If the vertices of the face do not belong all to the same boundary patch,
268 // the face is internal and is stored.
269 // Since each such face belongs to two elements, in order to ensure the
270 // face is stored only once, the vector vrV is converted into a set (ordered) used
271 // as the key value of the map, whose occurrence is thus unique.
272 if (! localcod) {
273 rankToNum(vrV);
274 set_n Face(vrV.begin(),vrV.end());
275 pair_nn nTnF(make_pair(itT->number(),i));
276 vector< pair_nn > VnTnF; VnTnF.push_back(nTnF);
277 retval = FI.insert(make_pair(Face, VnTnF));
278 if (! retval.second) {// add second element sharing this face
279 (retval.first->second).push_back(nTnF);
280 }
281 }
282 }
283 }
284 return FI;
285 }
286
287 //-------------------------------------------------------------------------------
288 // I/O utilities
289 //-------------------------------------------------------------------------------
290 /*!
291 Prints, on stream os, the definition of a mesh with vertices coordinates.
292 */
293 template<class T_>
printall(ostream & os) const294 void GeomFigureMesh<T_>::printall(ostream& os) const {
295 printInfo(os);
296
297 typename vector<T_>::const_iterator itT;
298 for (itT=listT_.begin(); itT != listT_.end(); itT++) {
299 os << "Element " << itT->number() << endl;
300 for (number_t j=1; j<=nb_main_vertices_by_element_; j++) {
301 listV_[itT->rankOfVertex(j)].print(os,TG_);
302 }
303 }
304 os << endl;
305
306 os << endl << "List of vertices :" << endl;
307 vector<Vertex>::const_iterator itV;
308 for (itV=listV_.begin(); itV != listV_.end(); itV++) {itV->print(os,TG_);}
309 os << endl;
310 }
311
312 /*!
313 Prints, on stream os, the definition of a mesh.
314 */
315 template<class T_>
print(ostream & os) const316 void GeomFigureMesh<T_>::print(ostream& os) const {
317 printInfo(os);
318
319 os << endl << "List of elements (vertices are given by their number):" << endl;
320 for (number_t nel=1; nel<=numberOfElements(); nel++) {
321 os << "Element " << nel << " : ";
322 vector<number_t> nV(element(nel));
323 for (size_t i=0; i<nV.size(); i++) {os << nV[i] << " ";}
324 os << endl;
325 }
326
327 os << endl << "List of elements (vertices are given by their rank in the vertex list):" << endl;
328 typename vector<T_>::const_iterator itT;
329 for (itT=listT_.begin(); itT != listT_.end(); itT++) {
330 itT->print(os);
331 }
332
333 os << endl << "List of vertices :" << endl;
334 vector<Vertex>::const_iterator itV;
335 number_t k;
336 for (itV=listV_.begin(), k=0; itV != listV_.end(); itV++,k++) {
337 os << "Rank " << setw(4) << k << " ";
338 itV->print(os,TG_);}
339 os << endl;
340 }
341
342 /*!
343 Prints, on stream ftex, Fig4TeX instructions to draw a mesh.
344 The observation direction is defined by the longitude Psi and the latitude Theta,
345 given in degrees. nbviews figures are drawn, all with the same latitude, the
346 longitude varying by a 360/nbviews increment.
347 */
348 template<class T_>
printTeX(ostream & ftex,const float Psi,const float Theta,const number_t nbviews,const string & DimProj,const bool withInterface,const bool withElems) const349 void GeomFigureMesh<T_>::printTeX(ostream& ftex, const float Psi, const float Theta,
350 const number_t nbviews, const string& DimProj,
351 const bool withInterface, const bool withElems) const {
352 ftex << "\\let\\showfigOne\\centerline" << endl;
353 ftex << "\\def\\showfigTwo#1#2{\\centerline{#1}\\nobreak\\medskip\\centerline{#2}}" << endl;
354 ftex << "\\input fig4tex.tex" << endl;
355 printTeXHeader(ftex);
356 ftex << "%" << endl;
357 ftex << "% 1. Definition of characteristic points" << endl;
358 ftex << "\\figinit{"<< DimProj << "}" << endl;
359
360 bool drawInterface = withInterface && numberOfInterfaces() > 0;
361 printTeXPoints(ftex,drawInterface);
362 ftex << "%" << endl;
363
364 float psi = Psi, deltapsi = 360./nbviews;
365 for (number_t numv=1; numv <= nbviews; numv++) {
366 ostringstream ss;
367 ss << "Subdiv. level "<< subdiv_level_
368 << ", long. "<< psi << "$^\\circ$, lat. " << Theta << "$^\\circ$";
369 if (drawInterface) {
370 createFileNBox(ftex,boundaryArea,psi,Theta,"A",string("Domain"));
371 createFileNBox(ftex,interfaceArea,psi,Theta,"B",string("Interfaces"));
372 ftex << "\\showfigTwo{\\box\\figBoxA\\hfil\\quad\\box\\figBoxB}{"
373 << ss.str() << "}" << endl;
374 }
375 else {
376 createFileNBox(ftex,boundaryArea,psi,Theta,"A",ss.str());
377 ftex << "\\showfigOne{\\box\\figBoxA}" << endl;
378 }
379 ftex << "%-------------------------------- End of figure --------------------------------" << endl;
380 ftex << "\\bigskip\\vfill" << endl;
381 psi += deltapsi;
382 }
383
384 // Draw the elements if requested (essentially usefull in 3D and to check high degree vertices)
385 if (withElems) {
386 ftex << "\\bigskip\\vfill\\eject" << endl;
387 ftex << "% Draw all the elements of the mesh" << endl;
388 printTeXInArea(ftex,0,subdomainArea);
389 ftex << "\\figdrawbegin{}"<< endl;
390 vector<number_t>::iterator itN;
391 number_t ind, indx(numberOfMainVerticesByElement());
392 // Draw all the elements of the mesh
393 number_t nbElts = numberOfElements();
394 for (number_t k=1; k<=nbElts; k++) {
395 vector<number_t> N = element(k);
396 ftex << "\\drawElem";
397 for (itN=N.begin(), ind=0; ind<indx; ind++,itN++) {
398 ftex << "{" << *itN << "}";
399 }
400 ftex << endl;
401 }
402 ftex << "\\figdrawend" << endl;
403 ftex << "\\figvisu{\\figBoxA}{" << nbElts << " elements of order " << order_ << "}{" << endl;
404 ftex << "% Write all the vertices as a whole" << endl;
405 ftex << "%\\figshowpts[1," << numberOfVertices() << "]" << endl;
406 ftex << "% Write all the vertices, element by element, including high order vertices if any" << endl;
407 ftex << "\\figset write(ptname={\\bf{#1}})" << endl;
408 // Write all the vertices, element by element, including high order vertices if any
409 for (number_t k=1; k<=nbElts; k++) {
410 vector<number_t> N = element(k);
411 itN=N.begin();
412 ftex << "\\figwritec[" << *itN++ ;
413 for (; itN<N.end(); itN++) {
414 ftex << "," << *itN ;
415 }
416 ftex << "]{}" << endl;
417 }
418 ftex << "}" << endl << "\\centerline{\\box\\figBoxA}" << endl;
419 }
420 ftex << "\\bye" << endl;
421 }
422
423 //-------------------------------------------------------------------------------
424 // Protected member functions
425 //-------------------------------------------------------------------------------
426 // defined here because of instantiation, but never used in practice
427 template<class T_>
algoSubdiv(const T_ & T,number_t & ElementNum,number_t & VertexNum,vector<T_> & listT,map_pair_num & SeenEdges)428 void GeomFigureMesh<T_>::algoSubdiv(const T_& T, number_t& ElementNum, number_t& VertexNum,
429 vector<T_>& listT, map_pair_num& SeenEdges) {
430 error("nofunc","GeomFigureMesh" ,"algoSubdiv");
431 }
432
433 //! build the mesh, compute high order vertices and check consistency
434 template<class T_>
buildNcheck(number_t & VertexNum)435 void GeomFigureMesh<T_>::buildNcheck(number_t& VertexNum) {
436 if (order_ < 1) {order_ = 1;}
437 // Subdivisions
438 buildMesh(VertexNum);
439 // Store the number of vertices of order 1
440 nb_main_vertices_ = listV_.size();
441
442 // End with computing high order vertices, if required, and some ckecks
443 try {
444 // Take order into account: compute additionnal high order vertices
445 if (order_ > 1) {
446 createHOV(VertexNum); // throw (HOV_error)
447 }
448 // Check number of vertices for each element and the number assigned to each of them
449 checkElements(); // throw (number_t *)
450 // Check the number assigned to each vertex
451 vector<Vertex>::const_iterator itV;
452 number_t k;
453 for (itV=listV_.begin(), k=minVertexNum_; itV != listV_.end(); itV++, k++) {
454 if (itV->number() != k) {
455 throw std::make_pair(itV->number(), k);
456 }
457 }
458 }
459 catch (HOV_error pb) {
460 cout << "*******************************************************" << endl;
461 cout << "From subdivision mesh constructor" << endl;
462 cout << " Error during high order vertices generation while processing" << endl;
463 cout << " element " << pb.elemNum << " and vertex " << pb.vertexNum << endl;
464 cout << pb.mes << endl;
465 cout << " Rank of vertex 1 = " << pb.rV1 << endl;
466 cout << " Rank of vertex 2 = " << pb.rV2 << endl;
467 cout << " Rank of vertex 3 = " << pb.rV3 << endl;
468 cout << "*******************************************************" << endl;
469 }
470 catch (number_t *values) {
471 cout << "*******************************************************" << endl;
472 cout << "From subdivision mesh constructor" << endl;
473 cout << " Error detected during post construction check:" << endl;
474 cout << "Computed element number : " << values[0] << endl;
475 cout << "Expected element number : " << values[1] << endl;
476 cout << values[2] << " vertices computed on this element " << endl;
477 cout << values[3] << " vertices expected." << endl;
478 cout << "*******************************************************" << endl;
479 }
480 catch (pair_nn values) {
481 cout << "*******************************************************" << endl;
482 cout << "From subdivision mesh constructor" << endl;
483 cout << " Error detected during post construction check:" << endl;
484 cout << "Computed vertex number : " << values.first << endl;
485 cout << "Expected vertex number : " << values.second << endl;
486 cout << "*******************************************************" << endl;
487 }
488 }
489
490 /*!
491 Build the mesh by successive subdivisions of the initial mesh.
492 Input arguments :
493 VertexNum : number of the last vertex created in the initial mesh
494 Output arguments :
495 VertexNum : number of the last vertex created in the mesh
496 */
497 template<class T_>
buildMesh(number_t & VertexNum)498 void GeomFigureMesh<T_>::buildMesh(number_t& VertexNum){
499 // At each subdivision, the elements are numbered starting from minElementNum_ (via
500 // the statement ElementNum=minElementNum_-1), since theses are all new elements.
501 // On the contrary, the vertices are all kept from the beginning and the vertex counter
502 // VertexNum must not be reinitialized.
503 // Nota: By construction, the elements stored in listT_ are grouped by subdomains.
504 for (number_t i=0; i<subdiv_level_; i++) {
505 number_t ElementNum = minElementNum_ - 1;
506 vector<T_> tmplistT;
507 // listT_.size() is the number of elements built during previous subdivision
508 tmplistT.reserve(subdivisionFactor_*listT_.size());
509 map_pair_num SeenEdges;
510 typename vector<T_>::iterator itT;
511 for (itT=listT_.begin(); itT != listT_.end(); itT++) {
512 algoSubdiv(*itT,ElementNum,VertexNum,tmplistT,SeenEdges);
513 }
514 listT_ = tmplistT;
515 }
516 initDefaultUserAttribute();
517 }
518 /*!
519 Check number of vertices for each element and the number assigned to each of them
520 */
521 template<class T_>
checkElements()522 void GeomFigureMesh<T_>::checkElements() {
523 number_t k;
524 typename vector<T_>::const_iterator itT;
525 for (itT=listT_.begin(), k=minElementNum_; itT != listT_.end(); itT++, k++) {
526 if (itT->numberOfVertices() != total_nb_vertices_by_element_ || itT->number() != k) {
527 number_t values[] = {itT->number(), k, itT->numberOfVertices(), total_nb_vertices_by_element_};
528 throw values;
529 }
530 }
531 }
532
533 /*!
534 Create high order vertices.
535 Input arguments :
536 \param VertexNum : number of the last vertex created so far in the mesh
537 Output arguments :
538 \param VertexNum : number of the last vertex created in the mesh by this procedure
539 */
540 template<class T_>
createHOV(number_t & VertexNum)541 void GeomFigureMesh<T_>::createHOV(number_t& VertexNum){
542 map_set_pair SeenEdges;
543 map_set_pair_in SeenFaces;
544 typename vector<T_>::iterator itT;
545
546 for (itT=listT_.begin(); itT != listT_.end(); itT++) {
547 // Vertices inside each edge
548 for (number_t numEdge=1; numEdge<=nb_edges_by_element_; numEdge++) {
549 createHOeV(*itT,order_,VertexNum,numEdge,SeenEdges);
550 }
551 // Vertices inside each face
552 for (number_t numFace=1; numFace<=nb_faces_by_element_; numFace++) {
553 createHOfV(*itT,order_,VertexNum,numFace,SeenFaces);
554 }
555 // Internal vertices
556 createHOiV(*itT,order_,VertexNum);
557 }
558 }
559 /*!
560 Create or retrieve high order vertices on the edge numEdge of element Elem.
561 Their numbers are consecutive on the edge.
562 The creation process of the new vertices and the numbering defined in
563 T_::numberingOfVertices are coherent.
564 */
565 template<class T_>
createHOeV(T_ & Elem,const number_t order,number_t & VertexNum,const number_t numEdge,map_set_pair & SeenEdges)566 void GeomFigureMesh<T_>::createHOeV(T_& Elem, const number_t order, number_t& VertexNum,
567 const number_t numEdge, map_set_pair& SeenEdges){
568 pair_nn prV = Elem.rkOfO1VeOnEdge(numEdge);
569 // Map key is a set in order to create a unique entry in the map, regardless of
570 // the way the edge is defined.
571 number_t rV1 = prV.first;
572 number_t rV2 = prV.second;
573 set_n Edge; Edge.insert(rV1); Edge.insert(rV2);
574 map_set_pair::iterator itSE;
575
576 if ((itSE=SeenEdges.find(Edge)) == SeenEdges.end()) {
577 // This edge has not yet been seen: compute the corresponding vertices
578 // Store the edge along with the corresponding couple (rank of the
579 // first new vertex created in the edge, rank of the first vertex
580 // defining the edge)
581 pair_nn NV (VertexNum, rV1);
582 SeenEdges.insert(make_pair(Edge,NV));
583
584 // Compute the new vertices along this edge
585 Vertex V1=listV_[rV1];
586 Vertex V2=listV_[rV2];
587 refnum_t localcod = V1.locCode() & V2.locCode();
588 // We compute the position of each point as a barycenter of the two end-points of the edge.
589 // The numbering is decreasing from V1 towards V2.
590 // V1 V2
591 // |----------|----------|----------|----------|----------|----------|----------|
592 // coef[0] = 1 2 . . . k-2 k-1
593 // coef[1] = k-1 k-2 . . . 2 1
594 vector<Point> VP(2);
595 VP[0] = V1.geomPt();
596 VP[1] = V2.geomPt();
597 real_t coef[2];
598 Point P;
599 for (number_t i1=1; i1<order; i1++) {
600 coef[0] = i1;
601 coef[1] = order-i1;
602 P = (this->*newVertexPt_)(localcod,coef,VP);
603 Elem.vertices_.push_back(VertexNum);// add rank of the new vertex
604 listV_.push_back(Vertex(++VertexNum,localcod,P));// store the new vertex
605 }
606 }
607 else {
608 // Edge found: get back the corresponding vertices whose numbers are consecutive
609 // along this edge, starting from num.
610 number_t num = (itSE->second).first;
611 // Find the edge orientation
612 if ( (itSE->second).second == rV1){
613 for (number_t i1=1; i1<order; i1++) {
614 Elem.vertices_.push_back(num++);// add rank of the vertex
615 }
616 }
617 else {
618 num += order-2;
619 for (number_t i1=1; i1<order; i1++) {
620 Elem.vertices_.push_back(num--);// add rank of the vertex
621 }
622 }
623 }
624 }
625 // defined here because of instantiation, but never used in practice
626 template<class T_>
printTeXHeader(ostream & ftex) const627 void GeomFigureMesh<T_>::printTeXHeader(ostream& ftex) const {
628 error("nofunc","GeomFigureMesh","printTeXHeader");
629 }
630 /*!
631 Prints, on stream ftex, Fig4TeX instructions to define points to be used for
632 the display of the mesh. These points are located on the boundaries of the domain
633 or on the interfaces.
634 */
635 template<class T_>
printTeXPoints(ostream & ftex,const bool withInterface) const636 void GeomFigureMesh<T_>::printTeXPoints(ostream& ftex, const bool withInterface) const {
637 printTeXInArea(ftex,0,boundaryArea);
638 if (withInterface) { printTeXInArea(ftex,0,interfaceArea); }
639 }
640
641 /*!
642 Prints, on stream ftex, Fig4TeX instructions to draw faces of a mesh belonging
643 to any area of kind TA. The faces are ordered from the farthest to the nearest
644 along the observation direction defined by the longitude psi and the latitude
645 theta, given in degrees.
646 */
647 template<class T_>
printTeXSortedAreaFoE(ostream & ftex,const topologicalArea TA,const float psi,const float theta) const648 void GeomFigureMesh<T_>::printTeXSortedAreaFoE(ostream& ftex, const topologicalArea TA,
649 const float psi, const float theta) const {
650 vector<TeXPolygon> BF;
651 number_t number_of_area = TG_.numberOf(TA);
652
653 // Collect all the boundary faces and keep their associated boundary number
654 for (number_t numArea=1; numArea<=number_of_area; numArea++) {
655 vector< vector<number_t> > FonB = rk_facesIn(TA,numArea);
656 vector< vector<number_t> >::iterator itFonB;
657 for (itFonB=FonB.begin(); itFonB != FonB.end(); itFonB++) {
658 BF.push_back(TeXPolygon(*itFonB,numArea,listV_));
659 }
660 // Define boundary colors
661 ftex << "\\def\\Color" << string(1,'@'+numArea) // ColorA, ColorB...
662 << "{" << colorOf(TA,numArea) << "}% " << TG_.nameOf(TA,numArea) << endl;
663 }
664 // Sort the faces from the farthest to the nearest
665 // 1. Compute the vector defining the observation direction
666 real_t psird = psi*pi_/180.;
667 real_t thetard = theta*pi_/180.;
668 real_t ct = std::cos(thetard);
669 // TeXPolygon::OD = Vect(ct*cos(psird), ct*sin(psird), sin(thetard));
670 TeXPolygon::initObsDir(Vect(ct*std::cos(psird), ct*std::sin(psird), std::sin(thetard)));
671 // 2. Apply sort algorithm
672 sort(BF.begin(),BF.end());
673 // 3. Print the sorted faces
674 if (number_of_area > 0) {
675 ftex << "% " << BF.size() << " faces on " << TG_.kindOf(TA) << " " << TG_.nameOf(TA,1);
676 }
677 for (number_t numArea=2; numArea<=number_of_area; numArea++) {
678 ftex << ", " << TG_.nameOf(TA,numArea);
679 }
680 ftex << endl;
681 vector<TeXPolygon>::iterator itBF;
682 vector<number_t>::const_iterator itF;
683 for (itBF=BF.begin(); itBF != BF.end(); itBF++) {
684 ftex << "\\drawFace";
685 for (itF=itBF->Vrank().begin(); itF != itBF->Vrank().end(); itF++) {
686 ftex << "{" << listV_[*itF].number() << "}";
687 }
688 ftex << "{\\Color" << string(1,'@'+itBF->attrib()) << "}";
689 ftex << endl;
690 }
691 }
692
693 /*!
694 Returns the list of edges of a mesh belonging to area num of kind TA
695 or belonging to any area of this kind if num=0.
696 Edges are defined by the rank of their vertices in the internal list of vertices.
697 */
698 template<class T_>
rk_edgesIn(const topologicalArea TA,const number_t num) const699 vector< pair_nn > GeomFigureMesh<T_>::rk_edgesIn(const topologicalArea TA, const number_t num) const {
700 vector< pair_nn > EonB;
701 refnum_t sig = lCodeOf(TA,num);
702
703 typename vector<T_>::const_iterator itT;
704 for (itT=listT_.begin(); itT != listT_.end(); itT++) {
705 for (number_t i=1; i<=nb_edges_by_element_; i++) {
706 pair_nn prV = itT->rkOfO1VeOnEdge(i);
707 refnum_t localcod = listV_[prV.first].locCode() & listV_[prV.second].locCode();
708 // If the edge belongs to the area, store it.
709 if (localcod & sig) {EonB.push_back(prV);}
710 }
711 }
712 return EonB;
713 }
714
715 /*!
716 Returns the list of faces of a mesh belonging to area num of kind TA
717 or belonging to any area of this kind if num=0.
718 Faces are defined by the rank of their vertices in the internal list of vertices.
719 */
720 template<class T_>
rk_facesIn(const topologicalArea TA,const number_t num) const721 vector< vector<number_t> > GeomFigureMesh<T_>::rk_facesIn(const topologicalArea TA, const number_t num) const {
722 vector< vector<number_t> > FonB;
723 refnum_t sig = lCodeOf(TA,num);
724
725 typename vector<T_>::const_iterator itT;
726 for (itT=listT_.begin(); itT != listT_.end(); itT++) {
727 for (number_t i=1; i<=nb_faces_by_element_; i++) {
728 vector<number_t> vrV = itT->rkOfO1VeOnFace(i);
729 vector<number_t>::const_iterator itvrV(vrV.begin());
730 refnum_t localcod = listV_[*itvrV].locCode();
731 while (++itvrV < vrV.end()) { localcod &= listV_[*itvrV].locCode(); }
732 // If the face belongs to the area, store it.
733 if (localcod & sig) {
734 if (itT->faceOrientation(i) > 0) { FonB.push_back(vrV); }
735 else {// orientation is negative : store vertices in the reverse order
736 FonB.push_back(vector<number_t> (vrV.rbegin(),vrV.rend())); }
737 }
738 }
739 }
740 return FonB;
741 }
742 //-------------------------------------------------------------------------------
743 // Private member functions
744 //-------------------------------------------------------------------------------
745 /*!
746 Write on stream ftex Fig4TeX instructions to create the graphical file and
747 the associated TeX box. Selected faces and edges belong to areas of kind TA.
748 Longitude psi and latitude theta define the observation direction.
749 boxName is the name of the TeX box created.
750 */
751 template<class T_>
createFileNBox(ostream & ftex,const topologicalArea TA,const float psi,const float theta,const char * boxName,const string & caption) const752 void GeomFigureMesh<T_>::createFileNBox(ostream& ftex, const topologicalArea TA,
753 const float psi, const float theta, const char* boxName,
754 const string& caption) const {
755 ftex << "\\figset proj(psi=" << psi << ", theta=" << theta << ")" << endl;
756 ftex << "% 2. Creation of the graphical file" << endl;
757 ftex << "\\figdrawbegin{}" << endl;
758 printTeXSortedAreaFoE(ftex,TA,psi,theta);
759 ftex << "\\figdrawend" << endl;
760 ftex << "%" << endl;
761 ftex << "% 3. Writing text on the figure" << endl;
762 ftex << "\\figvisu{\\figBox" << boxName << "}{" << caption << "}{" << endl;
763 // Display first points:
764 ftex << "\\figshowpts[1," << nb_main_vertices_by_element_ << "]" << endl;
765 ftex << "}" << endl;
766 }
767
768 // class instantiations:
769 template class GeomFigureMesh<Hexahedron>;
770 template class GeomFigureMesh<Tetrahedron>;
771 template class GeomFigureMesh<Quadrangle>;
772 template class GeomFigureMesh<Triangle>;
773
774 } // end of namespace subdivision
775 } // end of namespace xlifepp
776