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 Mesh.cpp
19 \authors D. Martin, Y. Lafranche, E. Lunéville, N. Kielbasiewicz
20 \since 26 may 2010
21 \date 5 mar 2014
22
23 \brief Implementation of xlifepp::Mesh class members and related functions
24 */
25
26 #include "Mesh.hpp"
27 #include "utils.h"
28
29 namespace xlifepp
30 {
31
operator <(const CrackData & cd) const32 bool CrackData::operator<(const CrackData& cd) const
33 {
34 if(dim != cd.dim) { return dim < cd.dim; }
35 else {return id < cd.id; }
36 }
37
computeNormal() const38 std::vector<real_t> CrackData::computeNormal() const
39 {
40 std::vector<real_t> normal;
41 if(dim == 1)
42 {
43 Point p=minimalBox.boundPt(2)-minimalBox.boundPt(1);
44 normal.resize(p.size());
45 if(p.x() == 0)
46 {
47 normal[0]=0.;
48 normal[1]=-p.z();
49 normal[2]=p.y();
50 return normal;
51 }
52 if(p.y() == 0)
53 {
54 normal[0]=p.z();
55 normal[1]=0.;
56 normal[2]=-p.x();
57 return normal;
58 }
59 if(p.z() == 0)
60 {
61 normal[0]=-p.y();
62 normal[1]=p.x();
63 normal[2]=0.;
64 return normal;
65 }
66 error("abnormal_failure");
67 return normal;
68 }
69 if(dim == 2)
70 {
71 Point p1=minimalBox.boundPt(2)-minimalBox.boundPt(1);
72 Point p2=minimalBox.boundPt(3)-minimalBox.boundPt(1);
73 return crossProduct(p1,p2);
74 }
75 error("nosuchcase","CrackData::computeNormal","dim=3");
76 return normal;
77 }
78
print(std::ostream & os) const79 void CrackData::print(std::ostream& os) const
80 {
81 if(theVerboseLevel <= 0) { return; }
82 os << "Cracked domain: id " << id << ", name " << domainName << ", dim " << dim << std::endl;
83 }
84
operator <<(std::ostream & os,const CrackData & crackData)85 std::ostream& operator<<(std::ostream& os, const CrackData& crackData)
86 {
87 crackData.print(os);
88 return os;
89 }
90
findId(std::vector<CrackData>::const_iterator it_b,std::vector<CrackData>::const_iterator it_e,number_t id)91 short int findId(std::vector<CrackData>::const_iterator it_b, std::vector<CrackData>::const_iterator it_e, number_t id)
92 {
93 std::vector<CrackData>::const_iterator it_cd;
94 number_t i=0;
95 for(it_cd=it_b; it_cd != it_e; ++it_cd, ++i)
96 {
97 if(it_cd->id == id) { return i; }
98 }
99 return -1;
100 }
101
102 //! void constructor
Mesh()103 Mesh::Mesh()
104 : geometry_p(0), lastIndex_(0), name_(""), comment_(""), isMadeOfSimplices_(true), order_(0), firstOrderMesh_p(0)
105 {}
106
107 //! construct mesh from file
Mesh(const string_t & filename,const string_t & na,IOFormat mft,number_t nodesDim)108 Mesh::Mesh(const string_t& filename, const string_t& na, IOFormat mft, number_t nodesDim)
109 : geometry_p(0), lastIndex_(0), name_(na), firstOrderMesh_p(0)
110 {
111
112 switch(mft)
113 {
114 case _msh:
115 loadGmsh(filename, nodesDim);
116 break;
117 case _mel:
118 loadMelina(filename, nodesDim);
119 break;
120 case _geo:
121 loadGeo(filename, nodesDim);
122 break;
123 case _ply:
124 loadPly(filename, nodesDim);
125 break;
126 case _vtk:
127 loadVtk(filename, nodesDim);
128 break;
129 case _vtu:
130 loadVtu(filename, nodesDim);
131 break;
132 case _medit:
133 loadMedit(filename, nodesDim);
134 break;
135 default:
136 error("not_yet_implemented", "Mesh::Mesh(const string_t&, const string_t&[, IOFormat[, Number]])");
137 break;
138 }
139 if(geometry_p==0) geometry_p = new Geometry(BoundingBox(computeBB()), "read from file " + filename, _fromFile);
140 }
141
copy(const Mesh & m)142 void Mesh::copy(const Mesh& m)
143 {
144 nodes = m.nodes;
145 copyAllButNodes(m);
146 }
147
copyAllButNodes(const Mesh & m)148 void Mesh::copyAllButNodes(const Mesh& m)
149 {
150 name_ = m.name_;
151 comment_ = m.comment_ + " (copy of " + name_ + ")";
152 geometry_p = m.geometry_p->clone();
153 vertices_ = m.vertices_;
154 order_ = m.order_;
155 isMadeOfSimplices_ = m.isMadeOfSimplices_;
156 lastIndex_ = m.lastIndex_;
157
158 // full copy of elements (pointers have to be updated)
159 elements_.resize(m.elements_.size());
160 std::vector<GeomElement*>::const_iterator ite;
161 std::vector<GeomElement*>::iterator it = elements_.begin();
162 number_t k = 0;
163 for(ite = m.elements_.begin(); ite != m.elements_.end(); ite++, it++, k++)
164 {
165 const GeomElement& elt = **ite;
166 GeomElement* nelt = new GeomElement(elt);
167 nelt->meshP() = this;
168 if(nelt->hasMeshElement()) // update MeshElement attributes (nodes, ...)
169 {
170 std::vector<Point*>::iterator itp;
171 std::vector<number_t>::iterator itn = nelt->meshElement()->nodeNumbers.begin();
172 for(itp = nelt->meshElement()->nodes.begin(); itp != nelt->meshElement()->nodes.end(); itp++, itn++)
173 *itp = &nodes[*itn-1];
174 nelt->meshElement()->geomMapData_p = 0;
175 }
176 else { error("geoelt_not_meshelement"); }
177 *it = nelt; //store new element
178 }
179
180 // full copy of mesh domains
181 domains_.resize(m.domains_.size());
182 for(number_t i = 0; i < m.domains_.size(); i++)
183 {
184 GeomDomain& dom = *m.domains_[i];
185 switch(dom.domType())
186 {
187 case _meshDomain:
188 {
189 GeomDomain* ndom = new GeomDomain(*this, dom.name(), dom.dim(), dom.description());
190 //copy elements
191 std::vector<GeomElement*>& elts = dom.meshDomain()->geomElements;
192 std::vector<GeomElement*>& nelts = ndom->meshDomain()->geomElements;
193 nelts.resize(elts.size());
194 std::vector<GeomElement*>::const_iterator ite;
195 std::vector<GeomElement*>::iterator itne = nelts.begin();
196 for(ite = elts.begin(); ite != elts.end(); ite++, itne++)
197 {
198 if((*ite)->isSideElement())
199 {
200 GeomElement* nelt = new GeomElement(**ite); //new side element
201 nelt->meshP()=this; //update parent mesh
202 std::vector<GeoNumPair>::iterator itg;
203 for(itg = nelt->parentSides().begin(); itg != nelt->parentSides().end(); itg++)
204 itg->first = elements_[itg->first->number()-1]; //update parent pointer
205 *itne = nelt;
206 }
207 else *itne = elements_[(*ite)->number()-1]; //copy of a mesh element
208 }
209
210 domains_[i] = ndom;
211 break;
212 }
213 case _compositeDomain:
214 case _analyticDomain:
215 default:
216 error("domain_notmesh",dom.name(), words("domain type",dom.domType()));
217 }
218 }
219
220 // optional vectors are not updated, reset to void
221 sides_.clear();
222 sideOfSides_.clear();
223 vertexElements_.clear();
224 firstOrderMesh_p = 0;
225 if(order_ == 1) firstOrderMesh_p = this;
226 setShapeTypes();
227 }
228
229 //! clear all pointers that have been copied at creation or copy
clear()230 void Mesh::clear()
231 {
232 if(geometry_p != 0) delete geometry_p;
233 geometry_p=0;
234 for(std::vector<GeomElement*>::iterator it = elements_.begin(); it != elements_.end(); it++)
235 if(*it != 0) delete *it;
236 for(std::vector<GeomElement*>::iterator it = sides_.begin(); it != sides_.end(); it++)
237 if(*it != 0) delete *it;
238 for(std::vector<GeomElement*>::iterator it = sideOfSides_.begin(); it != sideOfSides_.end(); it++)
239 if(*it != 0) delete *it;
240 if(firstOrderMesh_p != 0 && firstOrderMesh_p != this) delete firstOrderMesh_p;
241 for(std::vector<GeomDomain*>::iterator itd = domains_.begin(); itd != domains_.end(); itd++)
242 if(*itd != 0) delete *itd;
243 }
244
245 //! assign operator
operator =(const Mesh & m)246 Mesh& Mesh::operator=(const Mesh& m)
247 {
248 if(&m != this)
249 {
250 clear();
251 copy(m);
252 }
253 return *this;
254 }
255
256 //! access to a domain by its name
domain(const string_t & na) const257 const GeomDomain& Mesh::domain(const string_t& na) const
258 {
259 string_t sna = trim(na);
260 for(number_t i = 0; i < domains_.size(); i++)
261 {
262 if(domains_[i]->name() == sna) { return *domains_[i]; }
263 }
264 error("mesh_failfinddomain", na);
265 return *domains_[0]; //fake return
266 }
267
268 //! access to a domain by its name
domain(const string_t & na)269 GeomDomain& Mesh::domain(const string_t& na)
270 {
271 string_t sna = trim(na);
272 for(number_t i = 0; i < domains_.size(); i++)
273 {
274 if(domains_[i]->name() == sna) { return *domains_[i]; }
275 }
276 error("mesh_failfinddomain", na);
277 return *domains_[0]; //fake return
278 }
279
280 //! access to a domain by its number
domain(number_t n) const281 const GeomDomain& Mesh::domain(number_t n) const
282 {
283 if(n >= domains_.size()) { error("mesh_finddomainout", n, domains_.size()); }
284 return *domains_[n];
285 }
286
287 //! access to a domain by its number
domain(number_t n)288 GeomDomain& Mesh::domain(number_t n)
289 {
290 if(n >= domains_.size()) { error("mesh_finddomainout", n, domains_.size()); }
291 return *domains_[n];
292 }
293
294 //! access to a domain by its name
domainNumber(const string_t & na) const295 number_t Mesh::domainNumber(const string_t& na) const
296 {
297 string_t sna = trim(na);
298 for(number_t i = 0; i < domains_.size(); i++)
299 {
300 if(domains_[i]->name() == sna) { return i; }
301 }
302 error("mesh_failfinddomain", na);
303 return 0; //fake return
304 }
305
306 //! access to a domain by its name
domainNumber(const string_t & na)307 number_t Mesh::domainNumber(const string_t& na)
308 {
309 string_t sna = trim(na);
310 for(number_t i = 0; i < domains_.size(); i++)
311 {
312 if(domains_[i]->name() == sna) { return i; }
313 }
314 error("mesh_failfinddomain", na);
315 return 0; //fake return
316 }
317
318 //! access to a domain by its number
domainName(number_t n) const319 string_t Mesh::domainName(number_t n) const
320 {
321 if(n >= domains_.size()) { error("mesh_finddomainout", n, domains_.size()); }
322 return domains_[n]->name();
323 }
324
325 //! access to a domain by its number
domainName(number_t n)326 string_t Mesh::domainName(number_t n)
327 {
328 if(n >= domains_.size()) { error("mesh_finddomainout", n, domains_.size()); }
329 return domains_[n]->name();
330 }
331
332 //! try to find a MeshGenerator for a geometry and an element shape
defaultMeshGenerator(const Geometry & g,ShapeType sh)333 MeshGenerator Mesh::defaultMeshGenerator(const Geometry& g, ShapeType sh)
334 {
335 switch(g.dim())
336 {
337 case 1:
338 switch(g.shape())
339 {
340 case _segment:
341 case _setofpoints:
342 return _structured;
343 default:
344 return _gmsh;
345 }
346 case 2:
347 switch(g.shape())
348 {
349 case _parallelogram:
350 case _rectangle:
351 case _square:
352 return _structured;
353 case _disk:
354 case _setofelems:
355 return _subdiv;
356 default:
357 return _gmsh;
358 }
359 case 3:
360 switch(g.shape())
361 {
362 case _parallelepiped:
363 case _cuboid:
364 case _cube:
365 if(sh == _hexahedron) { return _structured; }
366 else { return _subdiv; }
367 case _ball:
368 case _revTrunk:
369 case _revCylinder:
370 case _revCone:
371 return _subdiv;
372 default:
373 return _gmsh;
374 }
375 }
376 return _gmsh;
377 }
378
379 /*---------------------------------------------------------------------------------------------
380 create test mesh with only one reference element (order one)
381 ---------------------------------------------------------------------------------------------*/
Mesh(ShapeType shape)382 Mesh::Mesh(ShapeType shape)
383 : name_("one reference "+words("shape",shape)+" element mesh")
384 {
385 number_t nbpt=0;
386 dimen_t dim=0;
387 BoundingBox bb;
388
389 switch(shape)
390 {
391 case _segment:
392 dim=1;
393 nbpt=2;
394 bb=BoundingBox(0,1);
395 nodes.resize(nbpt);
396 nodes[0] = Point(1.);
397 nodes[1] = Point(0.);
398 geometry_p = new Segment(_v1=nodes[0],_v2=nodes[1]);
399 break;
400
401 case _triangle:
402 dim=2;
403 nbpt=3;
404 bb=BoundingBox(0,1,0,1);
405 nodes.resize(nbpt);
406 nodes[0] = Point(1.,0.);
407 nodes[1] = Point(0.,1.);
408 nodes[2] = Point(0.,0.);
409 geometry_p = new Triangle(_v1=nodes[0],_v2=nodes[1],_v3=nodes[2]);
410 break;
411
412 case _quadrangle:
413 dim=2;
414 nbpt=4;
415 bb=BoundingBox(0,1,0,1);
416 nodes.resize(nbpt);
417 nodes[0] = Point(1.,0.);
418 nodes[1] = Point(1.,1.);
419 nodes[2] = Point(0.,1.);
420 nodes[3] = Point(0.,0.);
421 geometry_p = new Rectangle(_v1=nodes[3],_v2=nodes[0],_v4=nodes[2]);
422 break;
423
424 case _tetrahedron:
425 dim=3;
426 nbpt=4;
427 bb=BoundingBox(0,1,0,1,0,1);
428 nodes.resize(nbpt);
429 nodes[0] = Point(1.,0.,0.);
430 nodes[1] = Point(0.,1.,0.);
431 nodes[2] = Point(0.,0.,1.);
432 nodes[3] = Point(0.,0.,0.);
433 geometry_p = new Tetrahedron(_v1=nodes[0],_v2=nodes[1],_v3=nodes[2], _v4=nodes[3]);
434 break;
435
436 case _hexahedron:
437 dim=3;
438 nbpt=8;
439 bb=BoundingBox(0,1,0,1,0,1);
440 nodes.resize(nbpt);
441 nodes[0] = Point(1.,0.,0.);
442 nodes[1] = Point(1.,1.,0.);
443 nodes[2] = Point(0.,1.,0.);
444 nodes[3] = Point(0.,0.,0.);
445 nodes[4] = Point(1.,0.,1.);
446 nodes[5] = Point(1.,1.,1.);
447 nodes[6] = Point(0.,1.,1.);
448 nodes[7] = Point(0.,0.,1.);
449 geometry_p = new Cube(_v1=nodes[3],_v2=nodes[0],_v4=nodes[2],_v5=nodes[7]);
450 break;
451
452 case _prism:
453 dim=3;
454 nbpt=4;
455 bb=BoundingBox(0,1,0,1,0,1);
456 nodes.resize(nbpt);
457 nodes[0] = Point(1.,0.,0.);
458 nodes[1] = Point(0.,1.,0.);
459 nodes[2] = Point(0.,0.,0.);
460 nodes[3] = Point(0.,0.,1.);
461 geometry_p = new Prism(_v1=nodes[0],_v2=nodes[1],_v3=nodes[2],_direction=nodes[3]);
462 break;
463
464 case _pyramid:
465 dim=3;
466 nbpt=5;
467 bb=BoundingBox(0,1,0,1,0,1);
468 nodes.resize(nbpt);
469 nodes[0] = Point(0.,0.,0.);
470 nodes[1] = Point(1.,0.,0.);
471 nodes[2] = Point(1.,1.,0.);
472 nodes[3] = Point(0.,1.,0.);
473 nodes[4] = Point(0.,0.,1.);
474 geometry_p = new Pyramid(_v1=nodes[4],_v2=nodes[0],_v3=nodes[1],_v4=nodes[2],_apex=nodes[3]);
475 break;
476
477 default:
478 where("Mesh::Mesh(ShapeType shape)");
479 error("shape_not_handled",words("shape",shape));
480 }
481
482 Interpolation* interp_p = findInterpolation(Lagrange, standard, 1, H1);
483 RefElement* ref_p = findRefElement(shape, interp_p);
484 elements_.resize(1);
485 elements_[0] = new GeomElement(this, ref_p, dim, 1);
486 MeshElement* melt = elements_[0]->meshElement();
487 for(number_t i=0; i<nbpt; i++) melt->nodeNumbers[i]=i+1;
488 melt->vertexNumbers = melt->nodeNumbers;
489 melt->setNodes(nodes);
490 vertices_.resize(nodes.size());
491 for(number_t i = 0; i < nodes.size(); i++) { vertices_[i] = i + 1; }
492 MeshDomain* meshdom_p = (new GeomDomain(*this, "Omega", 3, "ref. "+words("shape",shape)))->meshDomain();
493 meshdom_p->geomElements = elements_;
494 meshdom_p->setShapeTypes();
495 domains_.push_back(meshdom_p);
496 lastIndex_ = 2;
497 buildGeomData();
498 setShapeTypes();
499 firstOrderMesh_p=this;
500 }
501
502 //! construction method from 1D geometry
build1DMesh(const Geometry & g,number_t order,MeshGenerator mg)503 void Mesh::build1DMesh(const Geometry& g, number_t order, MeshGenerator mg)
504 {
505 trace_p->push("Mesh::build1DMesh(const Geometry&[, Number[, MeshGenerator[, String]]]);");
506 if(g.dim() != 1) { error("geometry_only_1D"); }
507
508 lastIndex_=0;
509 firstOrderMesh_p=0;
510 geometry_p=g.clone();
511
512 if(mg == _defaultGenerator) { mg = defaultMeshGenerator(g); }
513
514 switch(mg)
515 {
516 case _gmsh:
517 {
518 string_t filename = "xlifepp_script.geo";
519 saveToGeo(*geometry_p, order, filename);
520 loadGeo(filename, geometry_p->dimPoint());
521 break;
522 }
523 case _structured:
524 if(order != 1) { error("bad_order", order, 1); }
525 switch(geometry_p->shape())
526 {
527 case _segment:
528 geometry_p->checkSideNamesAndUpdate(2);
529 {
530 Segment& geom(*geometry_p->segment());
531 meshP1Segment(geom, geom.n(), geom.sideNames());
532 }
533 break;
534 case _setofpoints:
535 geometry_p->checkSideNamesAndUpdate(2);
536 {
537 SetOfPoints& geom(*geometry_p->setofpoints());
538 meshP1Line(geom);
539 }
540 break;
541 default:
542 error("mesh_not_handled", geometry_p->domName(), words("shape",g.shape()));
543 break;
544 }
545 break;
546 default:
547 error("generator_not_handled", words("mesh generator", mg));
548 break;
549 }
550 trace_p->pop();
551 }
552
553 //! constructor from 2D and 3D geometries
buildMesh(const Geometry & g,ShapeType sh,number_t order,MeshGenerator mg,std::set<MeshOption> mos)554 void Mesh::buildMesh(const Geometry& g, ShapeType sh, number_t order, MeshGenerator mg, std::set<MeshOption> mos)
555 {
556 trace_p->push("Mesh::Mesh(const Geometry&, ShapeType[, Number[, MeshGenerator[, String]]]);");
557 if(g.dim() < 2)
558 {
559 warning("unused_shape_for_1D");
560 build1DMesh(g, order, mg);
561 }
562 else
563 {
564 std::vector<string_t> sideNames;
565 lastIndex_=0;
566 firstOrderMesh_p=0;
567 geometry_p=g.clone();
568
569 if(mg == _defaultGenerator) { mg = defaultMeshGenerator(g,sh); }
570
571 switch(mg)
572 {
573 case _gmsh:
574 {
575 string_t filename="xlifepp_script.geo";
576 saveToGeo(*geometry_p, sh, order, mos, filename);
577 loadGeo(filename, nodesDim(*geometry_p));
578 break;
579 }
580 case _structured:
581 if(order != 1) { error("bad_order",order,1); }
582 switch(geometry_p->shape())
583 {
584 case _parallelogram:
585 {
586 geometry_p->checkSideNamesAndUpdate(4);
587 Parallelogram& para= *geometry_p->parallelogram();
588 number_t n1 = para.n1(), n2=para.n2();
589 if(para.h().size()==4) {n1=number_t(para.length1()/para.h(1)); n2=number_t(para.length2()/para.h(1));}
590 switch(sh)
591 {
592 case _triangle : meshP1Parallelogram(para,n1,n2,para.sideNames()); break;
593 case _quadrangle: meshQ1Parallelogram(para,n1,n2,para.sideNames()); break;
594 default:
595 error("mesh_not_handled", geometry_p->domName(), words("shape",sh));
596 break;
597 }
598 }
599 break;
600
601 case _rectangle:
602 {
603 geometry_p->checkSideNamesAndUpdate(4);
604 Rectangle& rect= *geometry_p->rectangle();
605 number_t nx = rect.nx(), ny=rect.ny();
606 if(rect.h().size()==4) {nx=number_t(rect.xlength()/rect.h(1)); ny=number_t(rect.ylength()/rect.h(2));}
607 switch(sh)
608 {
609 case _triangle : meshP1Parallelogram(rect,nx,ny,rect.sideNames()); break;
610 case _quadrangle: meshQ1Parallelogram(rect,nx,ny,rect.sideNames()); break;
611 default:
612 error("mesh_not_handled", geometry_p->domName(), words("shape",sh));
613 break;
614 }
615 }
616 break;
617
618 case _square:
619 {
620 geometry_p->checkSideNamesAndUpdate(4);
621 Square& squa= *geometry_p->square();
622 number_t nx = squa.nx(), ny=squa.ny();
623 if(squa.h().size()==4) {nx=number_t(squa.xlength()/squa.h(1)); ny=number_t(squa.xlength()/squa.h(2));}
624 switch(sh)
625 {
626 case _triangle : meshP1Parallelogram(squa,nx,ny,squa.sideNames()); break;
627 case _quadrangle: meshQ1Parallelogram(squa,nx,ny,squa.sideNames()); break;
628 default:
629 error("mesh_not_handled", geometry_p->domName(), words("shape",sh));
630 break;
631 }
632 }
633 break;
634
635 case _parallelepiped:
636 geometry_p->checkSideNamesAndUpdate(6);
637 switch(sh)
638 {
639 case _tetrahedron:
640 meshP1Parallelepiped((*geometry_p->parallelepiped()), geometry_p->parallelepiped()->n1(),
641 geometry_p->parallelepiped()->n2(), geometry_p->parallelepiped()->n3(),
642 geometry_p->parallelepiped()->sideNames());
643 break;
644 case _hexahedron:
645 meshQ1Parallelepiped((*geometry_p->parallelepiped()), geometry_p->parallelepiped()->n1(),
646 geometry_p->parallelepiped()->n2(), geometry_p->parallelepiped()->n3(),
647 geometry_p->parallelepiped()->sideNames());
648 break;
649 case _prism:
650 meshPr1Parallelepiped((*geometry_p->parallelepiped()), geometry_p->parallelepiped()->n1(),
651 geometry_p->parallelepiped()->n2(), geometry_p->parallelepiped()->n3(),
652 geometry_p->parallelepiped()->sideNames());
653 break;
654 case _pyramid:
655 meshPy1Parallelepiped((*geometry_p->parallelepiped()), geometry_p->parallelepiped()->n1(),
656 geometry_p->parallelepiped()->n2(), geometry_p->parallelepiped()->n3(),
657 geometry_p->parallelepiped()->sideNames());
658 break;
659 default:
660 error("mesh_not_handled", geometry_p->domName(), words("shape",sh));
661 break;
662 }
663 break;
664
665 case _cuboid:
666 geometry_p->checkSideNamesAndUpdate(6);
667 switch(sh)
668 {
669 case _tetrahedron:
670 meshP1Parallelepiped((*geometry_p->cuboid()), geometry_p->cuboid()->nx(),
671 geometry_p->cuboid()->ny(), geometry_p->cuboid()->nz(), geometry_p->cuboid()->sideNames());
672 break;
673 case _hexahedron:
674 meshQ1Parallelepiped((*geometry_p->cuboid()), geometry_p->cuboid()->nx(),
675 geometry_p->cuboid()->ny(), geometry_p->cuboid()->nz(), geometry_p->cuboid()->sideNames());
676 break;
677 default:
678 error("mesh_not_handled", geometry_p->domName(), words("shape",sh));
679 break;
680 }
681 break;
682
683 case _cube:
684 geometry_p->checkSideNamesAndUpdate(6);
685 switch(sh)
686 {
687 case _tetrahedron:
688 meshP1Parallelepiped((*geometry_p->cube()), geometry_p->cube()->nx(),
689 geometry_p->cube()->ny(), geometry_p->cube()->nz(), geometry_p->cube()->sideNames());
690 break;
691 case _hexahedron:
692 meshQ1Parallelepiped((*geometry_p->cube()), geometry_p->cube()->nx(),
693 geometry_p->cube()->ny(), geometry_p->cube()->nz(), geometry_p->cube()->sideNames());
694 break;
695 default:
696 error("mesh_not_handled", geometry_p->domName(), words("shape", sh));
697 break;
698 }
699 break;
700
701 default:
702 error("mesh_not_handled", geometry_p->domName(), words("shape", sh));
703 break;
704 }
705 break;
706
707 case _subdiv:
708 switch(g.shape())
709 {
710 case _ball:
711 {
712 Ball& geom(*geometry_p->ball());
713 subdvMesh(geom, sh, geom.nbOctants(), geom.nbSubdiv(), order, geom.type(), geom.teXFilename());
714 break;
715 }
716 case _cube:
717 {
718 Cube& geom(*geometry_p->cube());
719 subdvMesh(geom, sh, geom.nbOctants(), geom.nbSubdiv(), order, geom.teXFilename());
720 break;
721 }
722 case _disk:
723 {
724 Disk& geom(*geometry_p->disk());
725 subdvMesh(geom, sh, geom.nbSubdiv(), order, geom.type(), geom.teXFilename());
726 break;
727 }
728 case _revCone:
729 {
730 RevCone& geom(*geometry_p->revCone());
731 subdvMesh(geom, sh, geom.nbSubdomains(), geom.nbSubdiv(), order, geom.type(), geom.teXFilename());
732 break;
733 }
734 case _revCylinder:
735 {
736 RevCylinder& geom(*geometry_p->revCylinder());
737 subdvMesh(geom, sh, geom.nbSubdomains(), geom.nbSubdiv(), order, geom.type(), geom.teXFilename());
738 break;
739 }
740 case _revTrunk:
741 {
742 RevTrunk& geom(*geometry_p->revTrunk());
743 subdvMesh(geom, sh, geom.nbSubdomains(), geom.nbSubdiv(), order, geom.type(), geom.teXFilename());
744 break;
745 }
746 case _setofelems:
747 {
748 SetOfElems& geom(*geometry_p->setofelems());
749 subdvMesh(geom.pts(), geom.elems(), geom.bounds(), geom.elemShape(), geom.nbSubdiv(), order, geom.teXFilename());
750 break;
751 }
752 default:
753 error("mesh_not_handled", geometry_p->domName(), words("shape", sh));
754 break;
755 }
756 break;
757
758 default:
759 error("generator_not_handled", words("mesh generator", mg));
760 break;
761 } // switch (mg)
762 } // else
763 trace_p->pop();
764 }
765
766 // constructor from a mesh msh to be subdivided nbsubdiv times
767 // msh is assumed to be of order 1 and made of elements all of the same shape.
768 // The resulting mesh is of order order.
Mesh(const Mesh & msh,number_t nbsubdiv,number_t order)769 Mesh::Mesh(const Mesh& msh, number_t nbsubdiv, number_t order)
770 :name_(msh.name_+"_subdivided")
771 {
772 using std::vector;
773 using std::set;
774 geometry_p=msh.geometry_p->clone();
775 if(msh.order_ > 1) { error("sub_mesh_order", msh.name_, msh.order_); }
776
777 // Fetch (unique) type of elements which must be the same for all the (non boundary) subdomains
778 ShapeType shty=_noShape;
779 for(vector<GeomDomain*>::const_iterator itdom=msh.domains_.begin(); itdom!=msh.domains_.end(); itdom++)
780 {
781 if((*itdom)->meshDomain()->isSideDomain()) { continue; }
782 else
783 {
784 const set<ShapeType>& sst=(*itdom)->meshDomain()->shapeTypes;
785 if(sst.size() != 1) { error("sub_mesh_elts", (*itdom)->meshDomain()->name(), msh.name_); }
786 if(shty != _noShape && shty != *(sst.begin())) { error("sub_mesh_shap", msh.name_); }
787 else { shty = *(sst.begin()); }
788 }
789 }
790 if(shty == _noShape) { error("sub_mesh_nodom", msh.name_); }
791
792 // Fetch boundaries
793 // Note: interface domains, if any, cannot be detected and thus are treated as boundaries, which would
794 // lead to an incorrect domain definition (in function build3Delements) since each interface
795 // element is shared by two elements of the main domain.
796 std::vector<std::vector<number_t> > bounds;
797 for(vector<GeomDomain*>::const_iterator itdom=msh.domains_.begin(); itdom!=msh.domains_.end(); itdom++)
798 {
799 if((*itdom)->meshDomain()->isSideDomain())
800 {
801 set<number_t> sn((*itdom)->meshDomain()->nodeNumbers());
802 bounds.push_back(vector<number_t>(sn.begin(), sn.end()));
803 }
804 }
805 // Fetch elements
806 short indHex[] = {3,7,2,6, 0,4,1,5}; // Vertices 4,8,3,7, 1,5,2,6 in XLiFE++ (cf. Hexahedron.cpp)
807 short indId[] = {0,1,2,3}; // Identity. OK for triangle, quadrangle and tetrahedron.
808 short* xl2sub = (shty==_hexahedron) ? indHex : indId;
809 std::vector<std::vector<number_t> > elems;
810 elems.reserve(elements_.size());
811 for(vector<GeomElement*>::const_iterator itel=msh.elements_.begin(); itel!=msh.elements_.end(); itel++)
812 {
813 vector<number_t> Vxl((*itel)->vertexNumbers()), Vsub(Vxl);
814 // Reordering of the vertices to match subdivision vertex numbering convention
815 short k=0;
816 for(vector<number_t>::iterator itVsub=Vsub.begin(); itVsub!=Vsub.end(); itVsub++)
817 {
818 *itVsub = Vxl[xl2sub[k++]];
819 }
820 elems.push_back(Vsub);
821 }
822 subdvMesh(msh.nodes, elems, bounds, shty, nbsubdiv, order);
823 // Retrieve original domain names. The order in which the domains are stored has not been changed.
824 /* No more true: The order in which the domains are stored may have changed.
825 Temporarily disabled. Everything is OK if the original domains have the same name as the new ones.
826 Otherwise, the user must check for the new names (with the printInfo() function) and uses them.
827 Fixing of this issue in progress.
828 for (vector<GeomDomain*>::const_iterator itdom=msh.domains_.begin(), itndom=domains_.begin(); itdom!=msh.domains_.end(); itdom++, itndom++)
829 {
830 (*itndom)->rename((*itdom)->name());
831 (*itndom)->setDescription((*itdom)->description());
832 }*/
833 }
834
835 /* -------------------------------------------------------------------------------------------
836 create a new mesh from a given one by spliting elements to other ones
837 can only split mesh having same elements
838 only few conversions are enable for the moment
839 hexahedron mesh -> pyramid mesh
840 --------------------------------------------------------------------------------------------*/
Mesh(const Mesh & mesh,ShapeType sh,const string_t name)841 Mesh::Mesh(const Mesh& mesh, ShapeType sh, const string_t name)
842 : name_(name)
843 {
844 trace_p->push("Mesh::Mesh(Mesh, ShapeType)");
845 std::set<ShapeType> shapeTypes; //!< list of element shape types in mesh domain (counted once)
846 std::vector<GeomElement*>::const_iterator itg;
847 for(itg = mesh.elements_.begin(); itg != mesh.elements_.end(); itg++) shapeTypes.insert((*itg)->shapeType());
848 if(shapeTypes.size()!=1) error("mesh_multiple_shapes");
849 ShapeType sh0=*shapeTypes.begin();
850 switch(sh0)
851 {
852 case hexahedron :
853 switch(sh)
854 {
855 case _pyramid : buildPyramidFromHexadron(mesh);
856 break;
857 default : error("mesh_bad_shape_conversion", words("shape",_hexahedron), words("shape",sh));
858 }
859 break;
860 default : error("mesh_bad_shape_conversion", words("shape",sh0), words("shape",sh));
861 }
862 trace_p->pop();
863 }
864
865 /*!
866 -------------------------------------------------------------------------------------------
867 Utility function used to complete a 1D mesh, after the nodes have been defined.
868 --------------------------------------------------------------------------------------------
869 */
complete1Dmesh(const string_t & domname,const std::vector<string_t> & sideNames)870 void Mesh::complete1Dmesh(const string_t& domname, const std::vector<string_t>& sideNames)
871 {
872 isMadeOfSimplices_ = true;
873 order_ = 1;
874 firstOrderMesh_p = this;
875
876 // construct mesh elements
877 Interpolation* interp_p = findInterpolation(Lagrange, standard, 1, H1);
878 RefElement* ref_p = findRefElement(_segment, interp_p);
879 dimen_t dim=nodes[0].size();
880 number_t nx = nodes.size();
881 elements_.resize(nx-1);
882 for(number_t k = 0; k < nx-1; k++)
883 {
884 elements_[k] = new GeomElement(this, ref_p, dim, k + 1);
885 MeshElement* melt = elements_[k]->meshElement();
886 melt->nodeNumbers[0] = k + 1;
887 melt->nodeNumbers[1] = k + 2;
888 melt->vertexNumbers = melt->nodeNumbers;
889 melt->setNodes(nodes);
890 }
891
892 // construct vertex indices (begins at 1)
893 vertices_.resize(nodes.size());
894 for(number_t i = 0; i < nodes.size(); i++) { vertices_[i] = i + 1; }
895
896 // construct main domain
897 string_t nmd=domname;
898 if(nmd=="") nmd="Omega";
899 MeshDomain* meshdom_p = (new GeomDomain(*this, nmd, 1, "segment [a,b]"))->meshDomain();
900 meshdom_p->geomElements = elements_;
901 domains_.push_back(meshdom_p);
902
903 // construct boundary domains
904 number_t k = nx;
905 if(sideNames[1] != "") //create domain {b}
906 {
907 meshdom_p = (new GeomDomain(*this, sideNames[1], 0, "point b of segment [a,b]"))->meshDomain();
908 meshdom_p->geomElements.push_back(new GeomElement(elements_[k - 2], 1, k));
909 k++;
910 domains_.push_back(meshdom_p);
911 }
912 if(sideNames[0] != "") //create domain {a}
913 {
914 meshdom_p = (new GeomDomain(*this, sideNames[0], 0, "point a of segment [a,b]"))->meshDomain();
915 meshdom_p->geomElements.push_back(new GeomElement(elements_[0], 2, k));
916 domains_.push_back(meshdom_p);
917 }
918 lastIndex_ = k;
919
920 //merge domains with same name
921 mergeDomainsWithSameName();
922
923 // compute measures and orientation of mesh elements
924 buildGeomData();
925 setShapeTypes();
926 }
927
928 /*---------------------------------------------------------------------------------------------
929 create regular mesh of segment [a,b]
930 sideNames[0]-> a, sideNames[1]-> b
931 if sideNames[i] is empty, the side is not a GeomDomain
932 nx is the number of nodes wanted on the segment (nx > 2)
933 ---------------------------------------------------------------------------------------------*/
meshP1Segment(Segment & seg,number_t nx,const std::vector<string_t> & sideNames)934 void Mesh::meshP1Segment(Segment& seg, number_t nx, const std::vector<string_t>& sideNames)
935 {
936 trace_p->push("Mesh::meshP1Segment");
937 Point a = seg.p1(), b = seg.p2();
938 dimen_t dim=a.size();
939 switch(dim)
940 {
941 case 3:
942 if(a(3)==0. && b(3)==0.)
943 {
944 if(a(2)==0. && b(2)==0.) { a.resize(1); b.resize(1); dim=1;} //Eric : force mesh in 1D
945 else { a.resize(2); b.resize(2); dim=2;} //Eric : force mesh in 2D
946 }
947 break;
948 case 2:
949 if(a(2)==0. && b(2)==0.) { a.resize(1); b.resize(1); dim=1;} //Eric : force mesh in 1D
950 break;
951 }
952 if(a == b) { error("void_elt", "segment"); }
953
954 // construct nodes
955 nodes.resize(nx);
956 for(number_t i = 0; i < nx; i++) { nodes[i] = a + i * (b-a)/(nx-1); }
957
958 complete1Dmesh(seg.domName(), sideNames);
959 trace_p->pop();
960 }
961
962 /*---------------------------------------------------------------------------------------------
963 Create a mesh of an opened polygonal line, made of consecutive segments whose end points are
964 two consecutive points taken from the input vector pts.
965 The dimension n of all the points should be the same and the mesh is created in Rn.
966 sideNames[0] is the domain name corresponding to the first point,
967 sideNames[1] is the domain name corresponding to the last point.
968 If sideNames[i] is empty, the side domain is not created.
969 ---------------------------------------------------------------------------------------------*/
meshP1Line(SetOfPoints & sp)970 void Mesh::meshP1Line(SetOfPoints& sp)
971 {
972 trace_p->push("Mesh::meshP1Line");
973 const std::vector<Point>& pts = sp.pts();
974 nodes.assign(pts.begin(), pts.end());
975 complete1Dmesh(sp.domName(),sp.sideNames());
976 trace_p->pop();
977 }
978
979 /*---------------------------------------------------------------------------------------------
980 create regular mesh of rectangle [a,b]x[c,d]
981 sideNames[0]->[a,b]x{c} sideNames[1]->{b}x[c,d] sideNames[2]->[a,b]x{d} sideNames[3]->{a}x[c,d]
982 if sideNames[i] is empty, the side is not a GeomDomain
983 nx and ny are the number of nodeswanted on each edge of the rectangle (nx > 2, ny > 2)
984 ---------------------------------------------------------------------------------------------*/
meshP1Parallelogram(Parallelogram & para,number_t nx,number_t ny,const std::vector<string_t> & sideNames)985 void Mesh::meshP1Parallelogram(Parallelogram& para, number_t nx, number_t ny,
986 const std::vector<string_t>& sideNames)
987 {
988 //initialisation
989 trace_p->push("Mesh::meshP1Parallelogram");
990 Point a = para.p(1), b = para.p(2), c = para.p(4);
991 dimen_t dim=a.size();
992 if(dim==3 && a(3) == 0. && b(3) == 0. && c(3) == 0.) //Eric : force 2D mesh
993 {
994 a.resize(2);
995 b.resize(2);
996 c.resize(2);
997 dim=2;
998 }
999 if(a == b || a == c) { error("degenerated_elt", "parallelogram"); }
1000
1001 isMadeOfSimplices_ = true;
1002 order_ = 1;
1003 firstOrderMesh_p = this;
1004
1005 //construct nodes
1006 nodes.resize(nx*ny);
1007 number_t k = 0;
1008 for(number_t j = 0; j < ny; j++)
1009 for(number_t i = 0; i < nx; i++, k++)
1010 nodes[k] = a+i*(b-a)/(nx-1)+j*(c-a)/(ny-1);
1011
1012 //construct mesh elements
1013 Interpolation* interp_p = findInterpolation(Lagrange, standard, 1, H1);
1014 RefElement* ref_p = findRefElement(_triangle, interp_p);
1015 elements_.resize(2 * (nx-1) * (ny-1));
1016 k = 0;
1017 for(number_t j = 0; j < ny-1 ; j++)
1018 {
1019 for(number_t i = 0; i < nx-1; i++)
1020 {
1021 elements_[k] = new GeomElement(this, ref_p, dim, k + 1);
1022 MeshElement* melt = elements_[k]->meshElement();
1023 number_t n = j * nx + i + 1;
1024 melt->nodeNumbers[0] = n + 1;
1025 melt->nodeNumbers[1] = n + nx;
1026 melt->nodeNumbers[2] = n;
1027 melt->vertexNumbers = melt->nodeNumbers; //update vertex numbers, same as node
1028 melt->setNodes(nodes); //update node pointers
1029 k++;
1030 elements_[k] = new GeomElement(this, ref_p, dim, k + 1);
1031 melt = elements_[k]->meshElement();
1032 melt->nodeNumbers[0] = n + nx;
1033 melt->nodeNumbers[1] = n + 1;
1034 melt->nodeNumbers[2] = n + nx + 1;
1035 melt->vertexNumbers = melt->nodeNumbers; //update vertex numbers
1036 melt->setNodes(nodes); //update node pointers
1037 k++;
1038 }
1039 }
1040
1041 //construct vertex indices (begins at 1)
1042 vertices_.resize(nodes.size());
1043 for(number_t i = 0; i < nodes.size(); i++) { vertices_[i] = i + 1; }
1044
1045 //construct main domain
1046 string_t nmd=para.domName();
1047 if(nmd=="") nmd="Omega";
1048 MeshDomain* meshdom_p = (new GeomDomain(*this, nmd, 2, "rectangle [a,b]x[c,d]"))->meshDomain();
1049 meshdom_p->geomElements = elements_;
1050 domains_.push_back(meshdom_p);
1051
1052 //construct boundary domains
1053 std::vector<string_t> descrip;
1054 descrip.push_back("boundary side [a,b]x{c} of rectangle [a,b]x[c,d]");
1055 descrip.push_back("boundary side {b}x[c,d] of rectangle [a,b]x[c,d]");
1056 descrip.push_back("boundary side [a,b]x{d} of rectangle [a,b]x[c,d]");
1057 descrip.push_back("boundary side {a}x[c,d] of rectangle [a,b]x[c,d]");
1058 for(number_t d = 0; d < sideNames.size(); d++)
1059 {
1060 if(sideNames[d] != "") //create a domain
1061 {
1062 meshdom_p = (new GeomDomain(*this, sideNames[d], 1, descrip[d]))->meshDomain();
1063 number_t nbelt = nx-1, numFace = 3;
1064 if(d == 1 || d == 3) {nbelt = ny-1; numFace = 2;}
1065 meshdom_p->geomElements.resize(nbelt);
1066 for(number_t e = 0; e < nbelt; e++)
1067 {
1068 number_t ne = 0;
1069 switch(d)
1070 {
1071 case 0:
1072 ne = 2 * e;
1073 break;
1074 case 1:
1075 ne = 2 * (e + 1) * (nx - 1) - 1;
1076 break;
1077 case 2:
1078 ne = 2 * (nx - 1) * (ny - 2) + 2 * e + 1;
1079 break;
1080 case 3:
1081 ne = 2 * e * (nx - 1);
1082 break;
1083 }
1084 meshdom_p->geomElements[e] = new GeomElement(elements_[ne], numFace, k + 1);
1085 k++;
1086 }
1087 domains_.push_back(meshdom_p);
1088 } // end of if (sideNames
1089 } // end of for (number_t d
1090 lastIndex_ = k ;
1091
1092 //merge domains with same name
1093 mergeDomainsWithSameName();
1094
1095 // compute measures and orientation of mesh elements
1096 buildGeomData();
1097 setShapeTypes();
1098
1099 // sides and sidesOfSides lists are not built (see buildSides and buildSideOfSides)
1100 trace_p->pop();
1101 }
1102
1103 //merge all domains declared in domains_ vector having the same name
mergeDomainsWithSameName()1104 void Mesh::mergeDomainsWithSameName()
1105 {
1106 if(domains_.size()==0) return; //no domain defined
1107 //list of domains by name
1108 std::map<string_t,std::vector<const GeomDomain*> > mapdom;
1109 std::map<string_t,std::vector<const GeomDomain*> >::iterator itm;
1110 std::vector<GeomDomain*>::iterator itd=domains_.begin();
1111 for(; itd!=domains_.end(); ++itd)
1112 {
1113 string_t dn=(*itd)->name();
1114 itm=mapdom.find(dn);
1115 if(itm==mapdom.end()) mapdom[dn]=std::vector<const GeomDomain*>(1,*itd);
1116 else itm->second.push_back(*itd);
1117 }
1118 //merge domains with same name
1119 domains_.clear();
1120 std::vector<const GeomDomain*>::iterator itc;
1121 for(itm=mapdom.begin(); itm!=mapdom.end(); ++itm)
1122 {
1123 if(itm->second.size()>1) //more than one, merge it
1124 {
1125 GeomDomain* gd = &xlifepp::merge(itm->second,itm->first);
1126 domains_.push_back(gd);
1127 for(itc=itm->second.begin(); itc!=itm->second.end(); ++itc) delete *itc; //delete old domains
1128 }
1129 else domains_.push_back(const_cast<GeomDomain*>(*itm->second.begin()));
1130 }
1131 }
1132
meshQ1Parallelogram(Parallelogram & para,number_t nx,number_t ny,const std::vector<string_t> & sideNames)1133 void Mesh::meshQ1Parallelogram(Parallelogram& para, number_t nx, number_t ny,
1134 const std::vector<string_t>& sideNames)
1135 {
1136 //initialisation
1137 trace_p->push("Mesh::meshQ1Parallelogram");
1138 Point a = para.p(1), b = para.p(2), c = para.p(4);
1139 dimen_t dim=a.size();
1140 if(dim==3 && a(3) == 0. && b(3) == 0. && c(3) == 0.) //Eric : force 2D mesh
1141 {
1142 a.resize(2);
1143 b.resize(2);
1144 c.resize(2);
1145 dim=2;
1146 }
1147 if(a == b || a == c) { error("degenerated_elt", "parallelogram"); }
1148
1149 isMadeOfSimplices_ = false;
1150 order_ = 1;
1151 firstOrderMesh_p = this;
1152
1153 //construct nodes
1154 nodes.resize(nx*ny);
1155 number_t k = 0;
1156 for(number_t j = 0; j < ny; j++)
1157 for(number_t i = 0; i < nx; i++, k++)
1158 nodes[k] = a+i*(b-a)/(nx-1)+j*(c-a)/(ny-1);
1159
1160 //construct mesh elements
1161 Interpolation* interp_p = findInterpolation(Lagrange, standard, 1, H1);
1162 RefElement* ref_p = findRefElement(_quadrangle, interp_p);
1163 elements_.resize((nx-1) * (ny-1));
1164 k = 0;
1165 for(number_t j = 0; j < ny-1 ; j++)
1166 {
1167 for(number_t i = 0; i < nx-1; i++, k++)
1168 {
1169 elements_[k] = new GeomElement(this, ref_p, dim, k + 1);
1170 MeshElement* melt = elements_[k]->meshElement();
1171 number_t n = j * nx + i + 1;
1172 melt->nodeNumbers[0] = n ;
1173 melt->nodeNumbers[1] = n + 1;
1174 melt->nodeNumbers[2] = n + nx + 1;
1175 melt->nodeNumbers[3] = n + nx;
1176 melt->vertexNumbers = melt->nodeNumbers; //update vertex numbers, same as node
1177 melt->setNodes(nodes); //update node pointers
1178 }
1179 }
1180
1181 //construct vertex indices (begins at 1)
1182 vertices_.resize(nodes.size());
1183 for(number_t i = 0; i < nodes.size(); i++) { vertices_[i] = i + 1; }
1184
1185 //construct main domain
1186 string_t nmd=para.domName();
1187 if(nmd=="") nmd="Omega";
1188 MeshDomain* meshdom_p = (new GeomDomain(*this, nmd, 2, "rectangle [a,b]x[c,d]"))->meshDomain();
1189 meshdom_p->geomElements = elements_;
1190 domains_.push_back(meshdom_p);
1191
1192 //construct boundary domains
1193 std::vector<string_t> descrip;
1194 descrip.push_back("boundary side [a,b]x{c} of rectangle [a,b]x[c,d]");
1195 descrip.push_back("boundary side {b}x[c,d] of rectangle [a,b]x[c,d]");
1196 descrip.push_back("boundary side [a,b]x{d} of rectangle [a,b]x[c,d]");
1197 descrip.push_back("boundary side {a}x[c,d] of rectangle [a,b]x[c,d]");
1198 for(number_t d = 0; d < sideNames.size(); d++)
1199 {
1200 if(sideNames[d] != "") //create a domain
1201 {
1202 meshdom_p = (new GeomDomain(*this, sideNames[d], 1, descrip[d]))->meshDomain();
1203 number_t nbelt = nx-1, numFace = 0;
1204 if(d == 1 || d == 3) nbelt = ny-1;
1205 switch(d)
1206 {
1207 case 0:
1208 numFace=1;
1209 break;
1210 case 1:
1211 numFace=2;
1212 break;
1213 case 2:
1214 numFace=3;
1215 break;
1216 case 3:
1217 numFace=4;
1218 break;
1219 }
1220 meshdom_p->geomElements.resize(nbelt);
1221 for(number_t e = 0; e < nbelt; e++)
1222 {
1223 number_t ne = 0;
1224 switch(d)
1225 {
1226 case 0:
1227 ne = e;
1228 break;
1229 case 1:
1230 ne = (e + 1) * (nx - 1) - 1;
1231 break;
1232 case 2:
1233 ne = (nx - 1) * (ny - 2) + e ;
1234 break;
1235 case 3:
1236 ne = e * (nx - 1);
1237 break;
1238 }
1239 meshdom_p->geomElements[e] = new GeomElement(elements_[ne], numFace, k + 1);
1240 k++;
1241 }
1242 domains_.push_back(meshdom_p);
1243 } // end of if (sideNames
1244 } // end of for (number_t d
1245 lastIndex_ = k ;
1246
1247 //merge domains with same name
1248 mergeDomainsWithSameName();
1249
1250 //compute measures and orientation of mesh elements
1251 buildGeomData();
1252 setShapeTypes();
1253
1254 //sides and sidesOfSides lists are not built (see buildSides and buildSideOfSides)
1255 trace_p->pop();
1256 }
1257
1258 /*---------------------------------------------------------------------------------------------
1259 create regular mesh of parallelepiped
1260 sideNames[0]-> [0,1]x[0,1]x{0} sideNames[1]-> [0,1]x[0,1]x{1}
1261 sideNames[2]-> [0,1]x{0}x[0,1] sideNames[3]-> [0,1]x{1}x[0,1]
1262 sideNames[4]-> {0}x[0,1]x[0,1] sideNames[5]-> {1}x[0,1]x[0,1]
1263 if sideNames[i] is empty, the side is not a identified as a GeomDomain
1264 nx, ny and nz are the number of intervals wanted on each edge of the parallelepiped (nx > 0, ny > 0, nz > 0)
1265 ---------------------------------------------------------------------------------------------*/
meshP1Parallelepiped(Parallelepiped & par,number_t nx,number_t ny,number_t nz,const std::vector<string_t> & sideNames)1266 void Mesh::meshP1Parallelepiped(Parallelepiped& par, number_t nx, number_t ny, number_t nz,
1267 const std::vector<string_t>& sideNames)
1268 { error("not_yet_implemented", "Mesh::meshP1Parallelepiped"); }
1269
meshQ1Parallelepiped(Parallelepiped & par,number_t nx,number_t ny,number_t nz,const std::vector<string_t> & sideNames)1270 void Mesh::meshQ1Parallelepiped(Parallelepiped& par, number_t nx, number_t ny, number_t nz,
1271 const std::vector<string_t>& sideNames)
1272 {
1273 // initialisation
1274 trace_p->push("Mesh::meshQ1Parallelepiped");
1275 Point a = par.p(1), b = par.p(2), c = par.p(4), d = par.p(5);
1276 if(a == b || a == c || a == d) { error("degenerated_elt", "parallelepiped"); }
1277
1278 isMadeOfSimplices_ = false;
1279 order_ = 1;
1280 firstOrderMesh_p = this;
1281
1282 // construct nodes
1283 nodes.resize(nx*ny*nz);
1284 number_t l = 0;
1285 for(number_t k = 0; k < nz; k++)
1286 for(number_t j = 0; j < ny; j++)
1287 for(number_t i = 0; i < nx; i++, l++)
1288 nodes[l] = a+i*(b-a)/(nx-1)+j*(c-a)/(ny-1)+k*(d-a)/(nz-1);
1289
1290 // construct mesh elements
1291 Interpolation* interp_p = findInterpolation(Lagrange, standard, 1, H1);
1292 RefElement* ref_p = findRefElement(_hexahedron, interp_p);
1293 elements_.resize((nx-1) * (ny-1) * (nz-1));
1294 l = 0;
1295 for(number_t k = 0; k < nz-1 ; k++)
1296 {
1297 for(number_t j = 0; j < ny-1 ; j++)
1298 {
1299 for(number_t i = 0; i < nx-1; i++, l++)
1300 {
1301 elements_[l] = new GeomElement(this, ref_p, 3, l + 1);
1302 MeshElement* melt = elements_[l]->meshElement();
1303 number_t q = nx * ny, n = k * q + j * nx + i + 1;
1304 melt->nodeNumbers[0] = n ;
1305 melt->nodeNumbers[1] = n + 1;
1306 melt->nodeNumbers[2] = n + nx + 1;
1307 melt->nodeNumbers[3] = n + nx;
1308 melt->nodeNumbers[4] = n + q;
1309 melt->nodeNumbers[5] = n + q + 1;
1310 melt->nodeNumbers[6] = n + nx + 1 + q;
1311 melt->nodeNumbers[7] = n + nx + q;
1312 melt->vertexNumbers = melt->nodeNumbers; // update vertex numbers, same as node
1313 melt->setNodes(nodes); // update node pointers
1314 }
1315 }
1316 }
1317
1318 // construct vertex indices (begins at 1)
1319 vertices_.resize(nodes.size());
1320 for(number_t i = 0; i < nodes.size(); i++) { vertices_[i] = i + 1; }
1321
1322 // construct main domain
1323 string_t nmd=par.domName();
1324 if(nmd=="") nmd="Omega";
1325 MeshDomain* meshdom_p = (new GeomDomain(*this, nmd, 3, "parallelepiped [a,b]x[c,d]x[e,f]"))->meshDomain();
1326 meshdom_p->geomElements = elements_;
1327 domains_.push_back(meshdom_p);
1328
1329 // construct boundary domains
1330 std::vector<string_t> descrip;
1331 descrip.push_back("boundary side [a,b]x[c,d]x{e} of parallelepiped [a,b]x[c,d]x[e,f]");
1332 descrip.push_back("boundary side [a,b]x[c,d]x{f} of parallelepiped [a,b]x[c,d]x[e,f]");
1333 descrip.push_back("boundary side [a,b]x{c}x[e,f] of parallelepiped [a,b]x[c,d]x[e,f]");
1334 descrip.push_back("boundary side [a,b]x{d}x[e,f] of parallelepiped [a,b]x[c,d]x[e,f]");
1335 descrip.push_back("boundary side {a}x[c,d]x[e,f] of parallelepiped [a,b]x[c,d]x[e,f]");
1336 descrip.push_back("boundary side {b}x[c,d]x[e,f] of parallelepiped [a,b]x[c,d]x[e,f]");
1337
1338 for(number_t d = 0; d < std::min(number_t(6), sideNames.size()); d++)
1339 {
1340 if(sideNames[d] != "") // create a domain
1341 {
1342 meshdom_p = (new GeomDomain(*this, sideNames[d], 2, descrip[d]))->meshDomain();
1343 number_t n1=0, n2=0, e = 0, numFace=0;
1344 switch(d)
1345 {
1346 case 0:
1347 n1=nx-1;
1348 n2=ny-1;
1349 numFace = 6;
1350 break;
1351 case 1:
1352 n1=nx-1;
1353 n2=ny-1;
1354 numFace = 3;
1355 break;
1356 case 2:
1357 n1=nx-1;
1358 n2=nz-1;
1359 numFace = 1;
1360 break;
1361 case 3:
1362 n1=nx-1;
1363 n2=nz-1;
1364 numFace = 4;
1365 break;
1366 case 4:
1367 n1=ny-1;
1368 n2=nz-1;
1369 numFace = 5;
1370 break;
1371 case 5:
1372 n1=ny-1;
1373 n2=nz-1;
1374 numFace = 2;
1375 break;
1376 }
1377 meshdom_p->geomElements.resize(n1 * n2);
1378 for(number_t j = 0; j < n2; j++)
1379 {
1380 for(number_t i = 0; i < n1; i++, e++)
1381 {
1382 number_t ne = 0;
1383 switch(d)
1384 {
1385 case 0:
1386 ne = e;
1387 break;
1388 case 1:
1389 ne = (nx - 1) * (ny - 1) * (nz - 2) + e ;
1390 break;
1391 case 2:
1392 ne = j * (nx - 1) * (ny - 1) + i;
1393 break;
1394 case 3:
1395 ne = (nx - 1) * (ny - 2) + j * (nx - 1) * (ny - 1) + i;
1396 break;
1397 case 4:
1398 ne = j * (nx - 1) * (ny - 1) + i * (nx - 1) ;
1399 break;
1400 case 5:
1401 ne = nx - 2 + j * (nx - 1) * (ny - 1) + i * (nx - 1) ;
1402 break;
1403 }
1404 meshdom_p->geomElements[e] = new GeomElement(elements_[ne], numFace, l + 1);
1405 l++;
1406 }
1407 }
1408 domains_.push_back(meshdom_p);
1409 } // end of if (sideNames
1410 } // end of for (number_t d
1411 lastIndex_ = l ;
1412
1413 //merge domains with same name
1414 mergeDomainsWithSameName();
1415
1416 // compute measures and orientation of mesh elements
1417 buildGeomData();
1418 setShapeTypes();
1419
1420 // sides and sidesOfSides lists are not built (see buildSides and buildSideOfSides)
1421 trace_p->pop();
1422 }
1423
1424 //mesh parallelepiped with prism, use extrusion of section
meshPr1Parallelepiped(Parallelepiped & par,number_t nx,number_t ny,number_t nz,const std::vector<string_t> & sideNames)1425 void Mesh::meshPr1Parallelepiped(Parallelepiped& par, number_t nx, number_t ny, number_t nz,
1426 const std::vector<string_t>& sideNames)
1427 {
1428 // initialisation
1429 trace_p->push("Mesh::meshPr1Parallelepiped");
1430 Point a = par.p(1), b = par.p(2), c = par.p(4), d = par.p(5);
1431 if(a == b || a == c || a == d) { error("degenerated_elt", "parallelepiped"); }
1432
1433 //mesh 2d section (a,c,d)
1434 Parallelogram S(_v1 = a, _v2 = c, _v4 = d, _nnodes = Numbers(nx,ny),
1435 _side_names=Strings(sideNames[0],sideNames[1],sideNames[2],sideNames[3]),_domain_name="Omega") ;
1436 Mesh section(S,_triangle,1,_structured);
1437
1438 //extrusion
1439 number_t ns=0, nb=0;
1440 if(sideNames[4]!="") ns=1;
1441 if(ns==0 && sideNames[5]!="") ns=1;
1442 for(number_t k=0; k<4 && nb==0; k++)
1443 if(sideNames[k]!="") nb=1;
1444 if(ns==0 && sideNames[5]!="") ns=1;
1445
1446 buildExtrusion(section, a, b, nz, 1, ns, nb);
1447 printInfo(thePrintStream);
1448
1449 //rename domains
1450 if(sideNames[4]!="") renameDomain("Omega_0",sideNames[4]);
1451 if(sideNames[5]!="") renameDomain("Omega_"+tostring(nz),sideNames[5]);
1452 string_t ndom=name_;
1453 if(ndom=="") ndom="Omega";
1454 renameDomain("Omega_e",ndom);
1455 std::set<string_t> sn;
1456 for(number_t k=0; k<4; k++) sn.insert(sideNames[k]);
1457 for(std::set<string_t>::iterator its =sn.begin(); its!=sn.end(); ++its)
1458 if(*its!="") renameDomain(*its+"_e",*its);
1459
1460 // compute measures and orientation of mesh elements
1461 buildGeomData();
1462 setShapeTypes();
1463
1464 printInfo(thePrintStream);
1465
1466 trace_p->pop();
1467 }
1468
1469 //mesh parallelepiped with pyramid
meshPy1Parallelepiped(Parallelepiped & par,number_t nx,number_t ny,number_t nz,const std::vector<string_t> & sideNames)1470 void Mesh::meshPy1Parallelepiped(Parallelepiped& par, number_t nx, number_t ny, number_t nz,
1471 const std::vector<string_t>& sideNames)
1472 {
1473 Mesh meshQ1(par, _hexahedron);
1474 buildPyramidFromHexadron(meshQ1);
1475 }
1476
1477 /*---------------------------------------------------------------------------------------------
1478 index all of side element of the mesh by travelling all sidedomains
1479 ---------------------------------------------------------------------------------------------*/
createSideEltIndex(std::map<string_t,GeomElement * > & sideEltIndex) const1480 void Mesh::createSideEltIndex(std::map<string_t, GeomElement*>& sideEltIndex) const //! index all the side elements of mesh
1481 {
1482 sideEltIndex.clear();
1483 for(number_t d = 0; d < nbOfDomains(); d++)
1484 {
1485 const MeshDomain* mdom = domains_[d]->meshDomain();
1486 if(mdom!=0 && mdom->isSideDomain()) xlifepp::createSideEltIndex(mdom->geomElements, sideEltIndex);
1487 }
1488 }
1489
1490 /*---------------------------------------------------------------------------------------------
1491 build the sides vector of all mesh elements (not built by default), side numbering not optimized
1492 side numberings of each mesh element are updated
1493 Be careful : this function has to be called after boundary domains construction
1494 ---------------------------------------------------------------------------------------------*/
buildSides()1495 void Mesh::buildSides()
1496 {
1497 trace_p->push("Mesh::buildSides");
1498 std::map<string_t, number_t> sideIndex; //temporary index where the string key is the ordered list of the side vertex numbers
1499 std::vector<GeomElement*>::iterator itel;
1500 std::vector<number_t>::iterator itn;
1501 std::map<string_t, number_t>::iterator itsn;
1502 number_t nbside = 0;
1503
1504 //first, initialize sideIndex map with sides from boundary domains
1505 for(number_t d = 1; d <= nbOfDomains(); d++)
1506 {
1507 const GeomDomain& dom = domains(d);
1508 if(dom.domType() == _meshDomain && dom.dim() == meshDim() - 1) //boundary domain of the mesh
1509 {
1510 std::vector<GeomElement*>::const_iterator itel;
1511 for(itel = dom.meshDomain()->geomElements.begin(); itel != dom.meshDomain()->geomElements.end(); itel++)
1512 {
1513 string_t key = (*itel)->encodeElement(); //encode vertex numbers
1514 if(sideIndex.find(key) == sideIndex.end()) //add GeomElement in map
1515 {
1516 sides_.push_back(*itel);
1517 nbside++;
1518 sideIndex[key] = nbside;
1519 }
1520 } // end of for (itel=
1521 } // end of if (dom
1522 } // end of for (number_t
1523
1524 // main loop on mesh elements
1525 for(itel = elements_.begin(); itel != elements_.end(); itel++)
1526 {
1527 MeshElement* melt = (*itel)->meshElement();
1528 for(number_t s = 0; s < (*itel)->numberOfSides(); s++) // loop on element sides
1529 {
1530 string_t key = (*itel)->encodeElement(s + 1); // encode vertex numbers of side s
1531 itsn = sideIndex.find(key); // search side in sideIndex map
1532 if(itsn == sideIndex.end()) // create a new side element
1533 {
1534 sides_.push_back(new GeomElement(*itel, s + 1, lastIndex_++)); //update sides_
1535 nbside++;
1536 sideIndex[key] = nbside; // update sideIndex
1537 melt->sideNumbers[s] = nbside; // update sideNumbers
1538 }
1539 else // update side numbering of element
1540 {
1541 melt->sideNumbers[s] = itsn->second; //update sideNumbers
1542 std::vector<GeoNumPair>& parside = sides_[itsn->second - 1]->parentSides(); //update parentSides_
1543 if(std::find(parside.begin(), parside.end(), GeoNumPair(*itel, s + 1)) == parside.end())
1544 { parside.push_back(GeoNumPair(*itel, s + 1)); }
1545 } // end of else
1546 } // end for (number_t s
1547 } //end for (itel
1548 trace_p->pop();
1549 }
1550
1551 /*---------------------------------------------------------------------------------------------
1552 build the sides of sides vector of all mesh elements (not built by default)
1553 note that side of side numbering of each mesh element is updated
1554 Important note : the parent element of a side of side (edge in 3D) is a side (face in 3D)
1555 thus the list of sides (sides_) has to be built before to built the list of sides of sides
1556 NOT CHECKED
1557 ---------------------------------------------------------------------------------------------*/
buildSideOfSides()1558 void Mesh::buildSideOfSides()
1559 {
1560 trace_p->push("Mesh::buildSideOfSides");
1561 if(meshDim() < 3) { error("mesh_nosos"); }
1562 if(sides_.size() == 0) { buildSides(); } //build list of sides if it is not built
1563
1564 std::map<string_t, number_t> sideIndex; //temporary index where the string key is the ordered list of the side of side vertex numbers
1565 std::vector<GeomElement*>::iterator itel;
1566 std::vector<number_t>::iterator itn;
1567 std::map<string_t, number_t>::iterator itsn;
1568 number_t nbside = 0;
1569
1570 // first, initialize sideIndex map with sides of sides from boundary domains of n-2 dimension
1571 for(number_t d = 1; d <= nbOfDomains(); d++)
1572 {
1573 const GeomDomain& dom = domains(d);
1574 if(dom.domType() == _meshDomain && dom.dim() == meshDim() - 2) //boundary domain of n-2 dimension
1575 {
1576 std::vector<GeomElement*>::const_iterator itel;
1577 for(itel = dom.meshDomain()->geomElements.begin(); itel != dom.meshDomain()->geomElements.end(); itel++)
1578 {
1579 string_t key = (*itel)->encodeElement(); //encode vertex numbers
1580 if(sideIndex.find(key) == sideIndex.end()) //add GeomElement in map
1581 {
1582 sides_.push_back(*itel);
1583 nbside++;
1584 sideIndex[key] = nbside;
1585 }
1586 } // end of for (itel=
1587 } // end of if (dom
1588 } // end of for (number_t
1589
1590 // main loop on sides element
1591 for(itel = sides_.begin(); itel != sides_.end(); itel++)
1592 {
1593 MeshElement* melt = (*itel)->parentSides()[0].first->meshElement(); // parent of side element as a mesh element
1594 for(number_t s = 0; s < (*itel)->geomRefElement()->nbSides(); s++) // loop on sides of side sit
1595 {
1596 string_t key = (*itel)->encodeElement(s); // encode vertex numbers of side s of side element itel
1597 itsn = sideIndex.find(key); // search side in sideIndex map
1598 if(itsn == sideIndex.end()) // create a new side element
1599 {
1600 sideOfSides_.push_back(new GeomElement(*itel, s + 1, lastIndex_++)); // update sideOfSides_
1601 nbside++;
1602 sideIndex[key] = nbside; // update sideIndex
1603 melt->sideOfSideNumbers.push_back(nbside); //update sideNumbers
1604 }
1605 else //update side numbering of element
1606 {
1607 melt->sideOfSideNumbers.push_back(itsn->second); // update sideOfSideNumbers
1608 sideOfSides_[itsn->second]->parentSides().push_back(GeoNumPair(*itel, s + 1)); // update parentSides_
1609 }
1610 } // end for (number_t s
1611 } // end for (itel
1612
1613 trace_p->pop();
1614 }
1615
1616 /*---------------------------------------------------------------------------------------------
1617 build for each vertex v of the mesh, the list of elements having v as vertex with its local number
1618 (not built by default)
1619 ---------------------------------------------------------------------------------------------*/
buildVertexElements()1620 void Mesh::buildVertexElements()
1621 {
1622 trace_p->push("Mesh::buildVertexElements");
1623 //build the list of global vertex number -> list of (elements,local vertex number)
1624 vertexElements_.resize(vertices_.size());
1625 std::vector<GeomElement*>::iterator itel; //mesh element iterator
1626 for(itel = elements_.begin(); itel != elements_.end(); itel++) //loop on mesh elements
1627 {
1628 for(number_t i = 1; i <= (*itel)->numberOfVertices(); i++)
1629 {
1630 vertexElements_[(*itel)->vertexNumber(i) - 1].push_back(GeoNumPair(*itel, i));
1631 }
1632 }
1633 trace_p->pop();
1634 }
1635
1636 /*! Update parent sides of elements involved in all side domains
1637 */
updateSideDomains()1638 void Mesh::updateSideDomains()
1639 {
1640 trace_p->push("Mesh::updateSideDomains");
1641 std::vector<GeomDomain*>::iterator itd=domains_.begin(); // iterator on list of geometric domains
1642 std::map<string_t,std::vector<GeoNumPair> > sideIndex; //temporary side index map
1643 for(; itd!=domains_.end(); ++itd) //loop on geometric domains
1644 {
1645 MeshDomain* mdom=(*itd)->meshDomain();
1646 if(mdom->isSideDomain()) //only side domains
1647 {
1648 if(sideIndex.size()==0) createSideIndex(sideIndex);
1649 mdom->updateSides(sideIndex);
1650 }
1651 }
1652 trace_p->pop();
1653 }
1654
1655 /*---------------------------------------------------------------------------------------------
1656 build geometric data
1657 ---------------------------------------------------------------------------------------------*/
buildGeomData()1658 void Mesh::buildGeomData()
1659 {
1660 if(geometry_p==0) geometry_p= new Geometry(BoundingBox(computeBB()), ""); //no geometry exists, create one to store BoundingBox
1661 else if(geometry_p->boundingBox.dim()==0) geometry_p->boundingBox=BoundingBox(computeBB()); // geometry exists, but void BoundingBox
1662 for(std::vector<GeomElement*>::iterator itElt = elements_.begin(); itElt != elements_.end(); itElt++)
1663 (*itElt)->meshElement()->computeMeasures();
1664 }
1665
1666 // set the shape types for all mesh domains
setShapeTypes()1667 void Mesh::setShapeTypes()
1668 {
1669 std::vector<GeomDomain*>::iterator itd=domains_.begin();
1670 for(; itd!=domains_.end(); itd++)
1671 {
1672 (*itd)->meshDomain()->setShapeTypes();
1673 }
1674 }
1675
1676
1677 /*---------------------------------------------------------------------------------------------
1678 mesh merging tool
1679 ---------------------------------------------------------------------------------------------*/
1680 /*!
1681 Merge a mesh to the current mesh with the following rules :
1682 - the two meshes must have the same dimensions (spaceDim and meshDim)
1683 - duplicated point are merged (duplicated means same points at a given tolerance, see Point class)
1684 - duplicated element (same order and same nodes) are merged
1685 notes : if meshes are not compatible, they are merged in a stupid way ...
1686 algorithm is not optimized !
1687 sides_, sideOfSides_, vertexElements_ vectors are not copied (recall buildSides)
1688
1689 mergeSharedBoundary : if true, same boundary domains are merged in a new one domain and original domains are kept
1690 with the naming rule : let na1 and na2 the first and second domain name
1691 the merged domain name is "na1 or na2"
1692 nmd : the optional name of the main domain (the whole domain), by default "main domain"
1693 */
merge(const Mesh & mesh,bool mergeSharedBoundary,const string_t & nmd)1694 void Mesh::merge(const Mesh& mesh, bool mergeSharedBoundary, const string_t& nmd)
1695 {
1696 trace_p->push("Mesh::merge");
1697 if(lastIndex_==0) //current mesh is presumed empty, do copy
1698 {
1699 copy(mesh);
1700 // renameDomain(0,nmd);
1701 return;
1702 }
1703
1704 if(mesh.spaceDim() != spaceDim()) { error("mesh_mismatch_dims", spaceDim(), mesh.spaceDim()); }
1705 if(mesh.meshDim() != meshDim()) { error("mesh_mismatch_dims", meshDim(), mesh.meshDim()); }
1706
1707 //update general data
1708 comment_ = "merge of " + name() + " and " + mesh.name();
1709 name_ += " + " + mesh.name();
1710 isMadeOfSimplices_ = isMadeOfSimplices_ && mesh.isMadeOfSimplices_;
1711 order_ = std::max(order_, mesh.order_);
1712 if(order_ > 2 && firstOrderMesh_p != 0) //destroy old underlying first order mesh
1713 {
1714 delete firstOrderMesh_p;
1715 firstOrderMesh_p = 0;
1716 }
1717 //clear sides, ... of current mesh (has to rebuilt later)
1718 sides_.clear();
1719 sideOfSides_.clear();
1720 vertexElements_.clear();
1721
1722 //update geometry attributes (bounding box)
1723 geometry_p->domName(geometry_p->domName() + " + " + mesh.geometry_p->domName());
1724 geometry_p->boundingBox += mesh.geometry_p->boundingBox;
1725
1726 real_t mesm=0.; //max of measure os elements
1727 std::vector<GeomElement*>::iterator it_elt1;
1728 for(it_elt1 = elements_.begin(); it_elt1 != elements_.end(); it_elt1++) //loop on elements of current mesh
1729 mesm=std::max(mesm,(*it_elt1)->meshElement()->measures[0]);
1730 real_t oldtol=Point::tolerance;
1731 if(mesm > 0) Point::tolerance =std::sqrt(mesm)/1000; //change tolerance in Point comparison
1732
1733 //merge nodes
1734 std::vector<Point>::iterator itnod1;
1735 std::vector<Point>::const_iterator itnod2;
1736 std::map<Point, number_t> newNodes; //new list of nodes with no duplicate node
1737 std::map<Point, number_t>::iterator itmpn; //iterator on previous map
1738 std::vector<number_t> mapNumbers(mesh.nodes.size(), 0); //numbering of nodes of mesh to add in new list of nodes
1739 number_t k = 1, k2 = 0;
1740 for(itnod1 = nodes.begin(); itnod1 != nodes.end(); itnod1++, k++) //insert nodes of current mesh
1741 newNodes.insert(std::pair<Point, number_t>(*itnod1, k));
1742 for(itnod2 = mesh.nodes.begin(); itnod2 != mesh.nodes.end(); itnod2++, k2++) //insert nodes of mesh to add
1743 {
1744 itmpn = newNodes.find(*itnod2);
1745 if(itmpn != newNodes.end()) // found
1746 mapNumbers[k2] = itmpn->second;
1747 else //not found, add in newNodes list
1748 {
1749 newNodes.insert(std::pair<Point, number_t>(*itnod2, k));
1750 mapNumbers[k2] = k;
1751 k++;
1752 }
1753 }
1754 Point::tolerance=oldtol; //reset old point tolerance
1755
1756 //update nodes of current mesh and clear newNodes map
1757 nodes.resize(newNodes.size());
1758 for(itmpn = newNodes.begin(); itmpn != newNodes.end(); itmpn++)
1759 { nodes[itmpn->second-1] = itmpn->first; }
1760 newNodes.clear();
1761
1762 //update vertices number (number are relative to nodes numbering)
1763 std::vector<number_t>::const_iterator itv;
1764 std::set<number_t> vert(vertices_.begin(),vertices_.end());
1765 for(itv = mesh.vertices_.begin(); itv != mesh.vertices_.end(); itv++) //loop on vertices of added mesh
1766 {
1767 number_t vn = mapNumbers[*itv-1];
1768 if(vert.find(vn)==vert.end()) vert.insert(vn); //vertex not in list
1769 }
1770 vertices_.assign(vert.begin(),vert.end());
1771 vert.clear();
1772
1773 //merge elements
1774 std::map<string_t, GeoNumPair> eltsIndex; //element index list of ordered vertex numbers -> (element,rank)
1775 std::map<number_t, number_t> eltNumbers; //numbering of elements of adding mesh in new list of elements
1776 std::vector<GeomElement*>::const_iterator it_elt2;
1777 std::map<string_t, GeoNumPair>::iterator it_msg; //iterator on eltsIndex map
1778
1779 k = 1;
1780 for(it_elt1 = elements_.begin(); it_elt1 != elements_.end(); it_elt1++, k++) //insert elements of current mesh
1781 {
1782 MeshElement* melt = (*it_elt1)->meshElement();
1783 melt->sideNumbers.clear(); //optionnal side informations are forgotten
1784 melt->sideOfSideNumbers.clear();
1785 for(number_t i = 0; i < melt->nodeNumbers.size(); i++) //update node pointers
1786 melt->nodes[i] = &(nodes[melt->nodeNumbers[i] - 1]);
1787 eltsIndex.insert(std::pair<string_t, GeoNumPair>((*it_elt1)->encodeElement(), GeoNumPair(*it_elt1, k)));
1788 }
1789 k2 = 1;
1790 for(it_elt2 = mesh.elements_.begin(); it_elt2 != mesh.elements_.end(); it_elt2++, k2++) //insert elements of mesh to add
1791 {
1792 //copy mesh element and update node numbers
1793 GeomElement* elt = new GeomElement(**it_elt2); //copy element
1794 elt->meshP() = this;
1795 MeshElement* melt = elt->meshElement();
1796 for(number_t i = 0; i < melt->nodeNumbers.size(); i++)
1797 melt->nodeNumbers[i] = mapNumbers[melt->nodeNumbers[i] - 1];
1798 for(number_t i = 0; i < melt->vertexNumbers.size(); i++)
1799 melt->vertexNumbers[i] = mapNumbers[melt->vertexNumbers[i] - 1];
1800
1801 string_t eltkey = elt->encodeElement();
1802 it_msg = eltsIndex.find(eltkey);
1803 if(it_msg != eltsIndex.end() && //same vertices
1804 it_msg->second.first->geomRefElement() == elt->geomRefElement()) //same reference element
1805 {
1806 eltNumbers[(*it_elt2)->number()] = it_msg->second.second;
1807 delete elt;
1808 }
1809 else //add element in current mesh
1810 {
1811 elt->number() = k; //update index of element
1812 melt->index() = k;
1813 for(number_t i = 0; i < melt->nodeNumbers.size(); i++) //update node pointers
1814 melt->nodes[i] = &(nodes[melt->nodeNumbers[i] - 1]);
1815 melt->sideNumbers.clear(); //optional side informations are not copied
1816 melt->sideOfSideNumbers.clear(); //too intricate, recall buildSides after
1817 elements_.push_back(elt);
1818 eltsIndex.insert(std::pair<string_t, GeoNumPair>(eltkey, GeoNumPair(elt, k)));
1819 eltNumbers[(*it_elt2)->number()] = k;
1820 k++;
1821 }
1822 }
1823
1824 //shift side element numbering of domains
1825 for(number_t d = 0; d < domains_.size(); d++)
1826 {
1827 MeshDomain* mdom = domains_[d]->meshDomain();
1828 for(it_elt2 = mdom->geomElements.begin(); it_elt2 != mdom->geomElements.end(); it_elt2++)
1829 {
1830 GeomElement* elt = *it_elt2;
1831 if(elt->isSideElement())
1832 {
1833 elt->number() = k;
1834 k++;
1835 }
1836 }
1837 }
1838
1839 number_t nbdom=domains_.size();
1840 //copy domains of adding mesh, domains may be renamed and pointers are updated but the domains are not merged !
1841 for(number_t d = 0; d < mesh.domains_.size(); d++)
1842 {
1843 MeshDomain* mdom = mesh.domains_[d]->meshDomain();
1844 //rename domain of adding mesh if its name already exists in current mesh
1845 string_t na = mdom->name();
1846 number_t n = 0;
1847 while(n < domains_.size())
1848 if(domains_[n]->name() != na) n++;
1849 else {na += "_2"; n = 0;}
1850 GeomDomain* newdom = new GeomDomain(*this, na, mdom->dim(), mdom->description());
1851 MeshDomain* ndom = newdom->meshDomain();
1852 for(it_elt2 = mdom->geomElements.begin(); it_elt2 != mdom->geomElements.end(); it_elt2++)
1853 {
1854 GeomElement* elt = *it_elt2;
1855 if(elt->isSideElement()) //create new side element
1856 {
1857 GeomElement* selt = new GeomElement(*elt);
1858 selt->number() = k;
1859 selt->meshP() = this;
1860 k++;
1861 selt->deleteMeshElement(); //destroy meshElement pointer to be safe
1862 for(number_t j = 0; j < selt->parentSides().size(); j++)
1863 {
1864 GeomElement* parelt = selt->parentSides()[j].first;
1865 if(parelt->isSideElement()) //is a side of side
1866 error("domain_type_not_handled", words("side of side"));
1867 else //is a side
1868 selt->parentSides()[j].first = elements_[eltNumbers[parelt->number()] - 1];
1869 }
1870 ndom->geomElements.push_back(selt);
1871 }
1872 else //not a side element
1873 ndom->geomElements.push_back(elements_[eltNumbers[elt->number()] - 1]);
1874 }
1875 domains_.push_back(ndom);
1876 }
1877
1878 dimen_t d=meshDim();
1879 number_t nbds = domains_.size();
1880 if(mergeSharedBoundary && nbds>nbdom) //merge shared boundary domains
1881 {
1882 //for each shared boundary domains pair (d1,d2), merge them in a new one with original name and rename the old one with n1 and n2
1883 for(number_t d1 = 0; d1 < nbdom; d1++) //loop on current mesh domains
1884 {
1885 MeshDomain* mdom1 = domains_[d1]->meshDomain();
1886 if(mdom1->dim() == d-1) //possible shared boundary
1887 {
1888 std::set<number_t> sn1=mdom1->nodeNumbers();
1889 for(number_t d2 = nbdom; d2 < nbds; d2++) //loop on added domains
1890 {
1891 MeshDomain* mdom2 = domains_[d2]->meshDomain();
1892 if(mdom2->dim() == d-1)
1893 {
1894 std::set<number_t> sn2=mdom2->nodeNumbers();
1895 if(sn1==sn2) //same domain
1896 {
1897 std::vector<GeomElement*>::const_iterator itg;
1898 std::map<string_t,GeomElement*> mgeo;
1899 for(itg = mdom1->geomElements.begin(); itg != mdom1->geomElements.end(); itg++)
1900 mgeo.insert(std::make_pair((*itg)->encodeElementNodes(),*itg));
1901 string_t rname = mdom1->name()+" or "+mdom2->name();
1902 GeomDomain* newdom = new GeomDomain(*this, rname, mdom1->dim(), mdom1->description());
1903 MeshDomain* ndom = newdom->meshDomain();
1904 for(itg = mdom2->geomElements.begin(); itg != mdom2->geomElements.end(); itg++)
1905 {
1906 std::map<string_t,GeomElement*>::iterator itm=mgeo.find((*itg)->encodeElementNodes());
1907 if(itm!=mgeo.end())
1908 {
1909 GeomElement* selt = new GeomElement(*itm->second);
1910 selt->number() = k;
1911 k++;
1912 selt->deleteMeshElement(); //destroy meshElement pointer to be safe
1913 selt->parentSides().push_back((*itg)->parentSides()[0]);
1914 ndom->geomElements.push_back(selt);
1915 }
1916 }
1917 if(ndom->geomElements.size()>0) domains_.push_back(ndom);
1918 else delete newdom; //void domain !!!
1919 }
1920 }
1921 }
1922 }
1923 }
1924 }
1925
1926 //construct main domain
1927 MeshDomain* meshdom_p = (new GeomDomain(*this, nmd, meshDim(), comment_))->meshDomain();
1928 //keep in main domain only elements of maximal dimension
1929 dimen_t maxdim=0;
1930 for(it_elt1 = elements_.begin(); it_elt1 != elements_.end(); it_elt1++)
1931 maxdim=std::max(maxdim, (*it_elt1)->meshElement()->elementDim());
1932 for(it_elt1 = elements_.begin(); it_elt1 != elements_.end(); it_elt1++)
1933 if((*it_elt1)->meshElement()->elementDim()==maxdim) meshdom_p->geomElements.push_back(*it_elt1);
1934 meshdom_p->setDescription(comment_+", made of "+tostring(meshdom_p->geomElements.size())+" elements");
1935 domains_.push_back(meshdom_p);
1936
1937 lastIndex_ = k;
1938 setShapeTypes();
1939 trace_p->pop();
1940 }
1941
1942 /*!
1943 create a underlying first order mesh (P1) if there exists some elements of order greater than one
1944 elements of order greater than one are split in elements of order one as usual :
1945 - P2 triangle split in 4 P1 triangles, P3 triangle split in 9 P1 triangles, ...
1946 - P2 tetrahedron split in 6 P1 tetrahedra, ...
1947 - Q2 hexahedron split in 6 P1 tetrahedra, ...
1948 */
createFirstOrderMesh() const1949 Mesh* Mesh::createFirstOrderMesh() const
1950 {
1951 Mesh* firstOrderMesh=0;
1952 if(order_==1 && isMadeOfSimplices_) {return firstOrderMesh;} //nothing is created
1953 trace_p->push("Mesh::createFirstOrderMesh"); //create new mesh
1954 firstOrderMesh=new Mesh;
1955 firstOrderMesh->name_=name_+"_first_order";
1956 firstOrderMesh->comment_="generated by Mesh::createFirstOrderMesh";
1957 firstOrderMesh->isMadeOfSimplices_=true;
1958 firstOrderMesh->order_=1;
1959 firstOrderMesh->lastIndex_=0;
1960 firstOrderMesh->nodes=nodes; //nodes are not modified when splitting mesh
1961 firstOrderMesh->geometry_p=geometry_p;
1962
1963 // build list of elements : keep first order element and split element of order greater than 1
1964 std::map<GeomElement*,std::vector<GeomElement*> > listel; //list of split elements
1965 typedef std::pair<GeomElement*,std::vector<GeomElement*> > mapair;
1966 std::vector<GeomElement*>::const_iterator itel;
1967 std::vector<GeomElement*>::iterator itel1;
1968 for(itel=elements_.begin(); itel!=elements_.end(); itel++)
1969 {
1970 GeomElement* elt=*itel;
1971 if(!elt->isSideElement()) //side element are not split
1972 {
1973 std::pair<std::map<GeomElement*,std::vector<GeomElement*> >::iterator,bool> ret;
1974 ret=listel.insert(mapair(elt,elt->splitP1())); //split element
1975 for(itel1=ret.first->second.begin(); itel1!=ret.first->second.end(); itel1++) //add to the list of elements of the first order mesh
1976 {
1977 firstOrderMesh->elements_.push_back(*itel1);
1978 }
1979 }
1980 }
1981 // update mesh pointer and number
1982 std::vector<GeomElement*>::iterator itel2;
1983 for(itel2=firstOrderMesh->elements_.begin(); itel2!=firstOrderMesh->elements_.end(); itel2++)
1984 {
1985 (*itel2)->meshP()=firstOrderMesh;
1986 (*itel2)->number()=firstOrderMesh->lastIndex_++;
1987 }
1988
1989 // update vertices (all nodes are vertices)
1990 number_t nbnodes=nodes.size();
1991 firstOrderMesh->vertices_.resize(nbnodes);
1992 std::vector<number_t>::iterator itv=firstOrderMesh->vertices_.begin();
1993 for(number_t i=1; i<=nbnodes; i++,itv++) { *itv=i; }
1994
1995 // update domains
1996 std::vector<GeomDomain*>::const_iterator itd;
1997 for(itd=domains_.begin(); itd!=domains_.end(); itd++)
1998 {
1999 GeomDomain* dom=*itd;
2000 if(dom->domType()==_meshDomain)
2001 {
2002 MeshDomain* meshdom = dom->meshDomain();
2003 MeshDomain* newdom = (new GeomDomain(*this, dom->name()+"_firstOrder", dom->dim(), dom->description()))->meshDomain();
2004 if(meshdom->isSideDomain()) //side domain, only side element
2005 {
2006 for(itel=meshdom->geomElements.begin(); itel!=meshdom->geomElements.end(); itel++)
2007 {
2008 GeomElement* elt=*itel;
2009 if(!elt->isSideElement())
2010 {
2011 where("Mesh::createFirstOrderMesh");
2012 error("side_elt_in_side_domain", meshdom->name(), elt->number());
2013 }
2014 std::vector<GeoNumPair>& parsides=elt->parentSides();
2015 std::vector<GeoNumPair>::iterator itgn;
2016 for(itgn=parsides.begin(); itgn!=parsides.end(); itgn++)
2017 {
2018 GeomElement* parelt=itgn->first;
2019 number_t side=itgn->second;
2020 if(parelt->isSideElement())
2021 {
2022 where("Mesh::createFirstOrderMesh");
2023 error("sos_elt_not_handled");
2024 }
2025 std::vector<std::pair<number_t,number_t> > sidenum=parelt->refElement()->splitP1Side(side);
2026 std::vector<std::pair<number_t,number_t> >::iterator its;
2027 for(its=sidenum.begin(); its!=sidenum.end(); its++)
2028 {
2029 GeomElement* eltP1=listel[parelt][its->first-1];
2030 GeomElement* sidelt=new GeomElement(eltP1,its->second,firstOrderMesh->lastIndex_++);
2031 newdom->geomElements.push_back(sidelt);
2032 }
2033 }
2034 }
2035 }
2036 else //plain domain, copy element
2037 {
2038 for(itel=meshdom->geomElements.begin(); itel!=meshdom->geomElements.end(); itel++)
2039 {
2040 for(itel1=listel[*itel].begin(); itel1!=listel[*itel].end(); itel1++)
2041 {
2042 newdom->geomElements.push_back(*itel1);
2043 }
2044 }
2045 }
2046 newdom->setShapeTypes();
2047 firstOrderMesh->domains_.push_back(newdom);
2048 }
2049 }
2050 trace_p->pop();
2051 return firstOrderMesh;
2052 }
2053
2054 /*--------------------------------------------------------------------------------------------
2055 remove a domain, be cautious when use it
2056 ---------------------------------------------------------------------------------------------*/
removeDomain(const string_t & na)2057 void Mesh::removeDomain(const string_t& na)
2058 {
2059 string_t sna = trim(na);
2060 number_t nd = domains_.size(), ina = nd;
2061 for(number_t i = 0; i < nd && ina==nd; i++)
2062 if(domains_[i]->name() == sna) ina=i;
2063 if(ina==nd) error("mesh_failfinddomain", na);
2064 delete domains_[ina];
2065 for(number_t i=ina; i<nd-1; i++) domains_[i]=domains_[i+1]; //shift
2066 domains_.resize(nd-1);
2067 }
2068
2069
2070 /*--------------------------------------------------------------------------------------------
2071 print mesh data
2072 ---------------------------------------------------------------------------------------------*/
printInfo(std::ostream & os) const2073 void Mesh::printInfo(std::ostream& os) const
2074 {
2075 os << "Mesh" << " '" << name_ << "'"; //name
2076 if(comment_ != "") { os << " (" << comment_ << ")"; } //comment
2077 os << "\n";
2078 os << " " << "space dimension" << " : " << spaceDim() << ", "
2079 << "element dimension" << " : " << meshDim() << "\n"; //global information
2080 os << " " << *geometry_p;
2081 os << " " << "number of elements" << " : " << nbOfElements() << ", "
2082 << "number of vertices" << " : " << nbOfVertices() << ", "
2083 << "number of nodes" << " : " << nbOfNodes() << ", "
2084 << "number of domains" << " : " << nbOfDomains();
2085 for(number_t i = 0; i < domains_.size(); i++)
2086 {
2087 os << "\n domain number " << i << ": " << domains_[i]->name() << " of dimension " << domains_[i]->dim()
2088 << ", made of " << domains_[i]->numberOfElements() << " elements (";
2089 if(domains_[i]->description().empty()) { os << "no description available"; }
2090 else { os << domains_[i]->description(); }
2091 os << ")";
2092 }
2093 os << "\n";
2094 }
2095
print(std::ostream & os) const2096 void Mesh::print(std::ostream& os) const
2097 {
2098 if(theVerboseLevel <= 0) { return; }
2099 printInfo(os);
2100
2101 if(theVerboseLevel < 3) { return; }
2102
2103 //list of elements, vertices, ... and nodes up to the verboselevel
2104 number_t s = elements_.size();
2105 os << "list of elements (" << s << ") :";
2106 number_t tvl = theVerboseLevel, n = std::min(tvl, s);
2107 verboseLevel(3);
2108 for(number_t i = 0; i < n; i++)
2109 { os << std::endl << (*elements_[i]); }
2110 if(n < s) { os << "\n...\n" << (*elements_[s - 1]) ; }
2111 verboseLevel(tvl);
2112
2113 s = vertices_.size();
2114 os << std::endl << words("list of") << " " << words("vertices") << " (" << s << ") :";
2115 n = std::min(tvl, s);
2116 for(number_t i = 0; i < n; i++)
2117 { os << "\n " << vertices_[i] << " -> " << nodes[vertices_[i] - 1]; }
2118 if(n < s) { os << "\n ...\n " << s << " -> " << nodes[vertices_[s - 1] - 1]; }
2119 verboseLevel(tvl);
2120
2121 s = nodes.size();
2122 os << std::endl << words("list of") << " " << words("nodes") << " (" << s << ") :";
2123 n = std::min(tvl, s);
2124 for(number_t i = 0; i < n; i++) { os << "\n " << i + 1 << " -> " << nodes[i]; }
2125 if(n < s) { os << "\n ...\n " << s << " -> " << nodes[s - 1]; }
2126 verboseLevel(tvl);
2127
2128 s = sides_.size();
2129 os << std::endl << words("list of") << " " << words("sides") << " (" << s << ") :";
2130 if(s == 0) { os << " unset"; }
2131 else
2132 {
2133 tvl = theVerboseLevel;
2134 n = std::min(tvl, s);
2135 verboseLevel(1);
2136 for(number_t i = 0; i < n; i++) { os << "\n side " << i + 1 << " -> " << (*sides_[i]); }
2137 if(n < s) { os << "\n ...\n " << "side " << s << " -> " << (*sides_[s - 1]); }
2138 verboseLevel(tvl);
2139 }
2140
2141 s = sideOfSides_.size();
2142 os << std::endl << words("list of") << " " << words("side of sides") << " (" << s << ") :";
2143 if(s == 0) { os << " unset"; }
2144 else
2145 {
2146 n = std::min(tvl, s);
2147 for(number_t i = 0; i < n; i++) { os << "\n " << words("side of side") << " " << i + 1 << " -> " << (*sideOfSides_[i]); }
2148 if(n < s) { os << "\n ...\n " << words("side of side") << " " << s << " -> " << (*sideOfSides_[s - 1]); }
2149 }
2150
2151 s = domains_.size();
2152 os << std::endl << words("list of") << " " << words("domains") << " (" << s << ") :";
2153 for(number_t i = 0; i < s; i++)
2154 { os << std::endl << (*domains_[i]); }
2155
2156 s = vertexElements_.size();
2157 os << std::endl << words("list of") << " " << words("elements") << " " << words("by vertex") << " (" << s << ") :";
2158 if(s == 0) { os << " unset"; }
2159 else
2160 {
2161 n = std::min(tvl, s);
2162 for(number_t i = 0; i < n; i++)
2163 {
2164 os << "\n " << words("vertex") << " " << vertices_[i] << " : ";
2165 for(number_t k = 0; k < vertexElements_[i].size(); k++)
2166 { os << " " << vertexElements_[i][k].first->number(); }
2167 }
2168 if(n < s)
2169 {
2170 os << "\n ...\n" << " " << words("vertex") << " " << vertices_[s - 1] << " : ";
2171 for(number_t k = 0; k < vertexElements_[s - 1].size(); k++)
2172 { os << " " << vertexElements_[s - 1][k].first->number(); }
2173 }
2174 }
2175 }
2176
operator <<(std::ostream & os,const Mesh & mesh)2177 std::ostream& operator<<(std::ostream& os, const Mesh& mesh)
2178 {
2179 mesh.print(os);
2180 return os;
2181 }
2182
addSuffix(const string_t & s)2183 void Mesh::addSuffix(const string_t& s)
2184 {
2185 name_+="_"+s;
2186 geometry_p->addSuffix(s);
2187 for(number_t i = 0; i < domains_.size(); ++i) { domains_[i]->addSuffix(s); }
2188 }
2189
2190 //! adaptative domains for levelset methods
adaptDomains(number_t mainDomNum,const std::vector<number_t> & elementsColorMap)2191 void Mesh::adaptDomains(number_t mainDomNum, const std::vector<number_t>& elementsColorMap)
2192 {
2193 // the domain number is the index in Mesh::domains_
2194 // nodesColorMap affects a domain number to each node of a mesh, even side domains
2195
2196 // 1. we detect which domains are plain domains
2197 // In the meantime, we build list of element numbers
2198 std::vector<number_t> plainDomains, sideDomains;
2199 std::map<number_t, std::set<number_t> > eltNumsPerDomains;
2200
2201 for(number_t d=0; d < domains_.size(); ++d)
2202 {
2203 if(domains_[d]->meshDomain()->isSideDomain()) { sideDomains.push_back(d); }
2204 else { plainDomains.push_back(d); }
2205 for(number_t e=0; e < domains_[d]->meshDomain()->numberOfElements(); ++e)
2206 {
2207 eltNumsPerDomains[d].insert(domains_[d]->meshDomain()->element(e+1)->number());
2208 }
2209 }
2210
2211 // 2. we build additional data between domains
2212
2213 // 2.a we build parent domains of side domains
2214 // it is enough to determine domains containing parents of the first element
2215 // so we have to update parents management
2216 updateSideDomains();
2217 std::map<number_t, std::set<number_t> > parentSideDomains;
2218 std::map<number_t, std::set<number_t> > boundaryPlainDomains;
2219
2220 for(number_t sd=0; sd < sideDomains.size(); ++sd)
2221 {
2222 number_t sdomId=sideDomains[sd];
2223 GeomElement* gelt=domains_[sdomId]->meshDomain()->element(1);
2224 std::vector<GeoNumPair>& parents=gelt->parentSides();
2225 std::vector<number_t> parentEltNums;
2226 for(number_t p=0; p < parents.size(); ++p) { parentEltNums.push_back(parents[p].first->number()); }
2227 for(number_t pe=0; pe < parentEltNums.size(); ++pe)
2228 {
2229 for(number_t pd=0; pd < plainDomains.size(); ++pd)
2230 {
2231 number_t pdomId=plainDomains[pd];
2232 if(eltNumsPerDomains[pdomId].count(parentEltNums[pe])) { parentSideDomains[sdomId].insert(pdomId); }
2233 }
2234 }
2235 }
2236
2237 // 2.b we build side domains of plain domains
2238 // we use previous data and reverse it
2239 std::map<number_t, std::set<number_t> >::const_iterator it_psd;
2240 for(it_psd=parentSideDomains.begin(); it_psd != parentSideDomains.end(); ++it_psd)
2241 {
2242 std::set<number_t>::const_iterator it;
2243 for(it=it_psd->second.begin(); it != it_psd->second.end(); ++it)
2244 {
2245 boundaryPlainDomains[*it].insert(it_psd->first);
2246 }
2247 }
2248
2249 // 3. we update plain domains by regenerating elements lists
2250 std::map<number_t,std::vector<GeomElement*> > eltsByDom;
2251
2252 // 3.a we first build list of elements per plain domain
2253 // as number of plain elements in mesh does not change,
2254 // plain elements numbering is left unchanged
2255 for(number_t e=0; e < elementsColorMap.size(); ++e)
2256 {
2257 eltsByDom[elementsColorMap[e]].push_back(elements_[e]);
2258 }
2259
2260 // 3.b we update plain domains
2261 for(number_t d=0; d < plainDomains.size(); ++d)
2262 {
2263 number_t domId=plainDomains[d];
2264 std::vector<GeomElement*>& v_ge =domains_[domId]->meshDomain()->geomElements;
2265 v_ge.resize(eltsByDom[domId].size());
2266 for(number_t e=0; e < eltsByDom[d].size(); ++e) { v_ge[e]=eltsByDom[domId][e]; }
2267 }
2268
2269 // 4 we update side domains by regenerating elements lists
2270
2271 // 4.a first, we get available numbering for side elements
2272 // To do so, we need the max index used and to keep numbers for plain elements and
2273 // for side elements that are not on interfaces between plain domains
2274 std::map<number_t, std::pair<number_t, number_t> > eltRangePerDomains;
2275 std::map<number_t, std::set<number_t> >::const_iterator it_enpd;
2276 for(it_enpd=eltNumsPerDomains.begin(); it_enpd != eltNumsPerDomains.end(); ++it_enpd)
2277 {
2278 eltRangePerDomains[it_enpd->first]=std::make_pair(*it_enpd->second.begin(), *it_enpd->second.rbegin());
2279 }
2280
2281 number_t maxNum=0;
2282 std::map<number_t, std::pair<number_t, number_t> >::const_iterator it_erpd;
2283 for(it_erpd=eltRangePerDomains.begin(); it_erpd != eltRangePerDomains.end(); ++it_erpd)
2284 {
2285 if(it_erpd->second.second > maxNum) { maxNum = it_erpd->second.second; }
2286 }
2287
2288 std::set<number_t> availableNumberingSet;
2289 for(number_t i=0; i < maxNum; ++i)
2290 {
2291 availableNumberingSet.insert(i+1);
2292 }
2293
2294 for(it_enpd=eltNumsPerDomains.begin(); it_enpd != eltNumsPerDomains.end(); ++it_enpd)
2295 {
2296 std::set<number_t>::const_iterator it_en;
2297 if(!(parentSideDomains.count(it_enpd->first) && parentSideDomains[it_enpd->first].size() > 1))
2298 {
2299 for(it_en=it_enpd->second.begin(); it_en != it_enpd->second.end(); ++it_en)
2300 {
2301 availableNumberingSet.erase(*it_en);
2302 }
2303 }
2304 }
2305
2306 std::vector<number_t> availableNumbering(availableNumberingSet.begin(), availableNumberingSet.end());
2307
2308 // 4.b then, we compute side domain of every internal plain domain
2309 // an internal plain domain is a domain that has only shared boundaries
2310 // in only one domain
2311 std::map<number_t, std::vector<GeoNumPair> > sideElts;
2312 for(number_t pd=0; pd < plainDomains.size(); ++pd)
2313 {
2314 number_t pdomNum=plainDomains[pd];
2315 if(pdomNum != mainDomNum)
2316 {
2317 // domain number of boundary
2318 // in this case boundaryPlainDomains[pdomNum].size() == 1
2319 number_t sDomNum=*boundaryPlainDomains[pdomNum].begin();
2320 std::map<string_t,std::vector<GeoNumPair> > sideIndex;
2321 std::map<string_t,std::vector<GeoNumPair> >::const_iterator it;
2322 xlifepp::createSideIndex(domains_[pdomNum]->meshDomain()->geomElements, sideIndex);
2323 for(it=sideIndex.begin(); it != sideIndex.end(); ++it)
2324 {
2325 if(it->second.size() == 1) // side element is not shared => it is a boundary element
2326 {
2327 sideElts[sDomNum].push_back(it->second[0]);
2328 }
2329 }
2330 }
2331 }
2332
2333 // 4.c we update side domains
2334 std::map<number_t, std::vector<GeoNumPair> >::const_iterator it_se;
2335 for(it_se=sideElts.begin(); it_se != sideElts.end(); ++it_se)
2336 {
2337 number_t sdomNum=it_se->first;
2338 // we have to clean list of elements of domain
2339 std::vector<GeomElement*>& v_ge =domains_[sdomNum]->meshDomain()->geomElements;
2340 for(std::vector<GeomElement*>::iterator it = v_ge.begin(); it != v_ge.end(); it++)
2341 if(*it != 0) delete *it;
2342
2343 // we build the new list by using available numbering
2344 v_ge.resize(it_se->second.size());
2345 for(number_t e=0; e < it_se->second.size(); ++e)
2346 {
2347 number_t eltNum;
2348 if(e < availableNumbering.size()) { eltNum=availableNumbering[e]; }
2349 else
2350 {
2351 eltNum=maxNum++;
2352 }
2353 v_ge[e]=new GeomElement(it_se->second[e].first, it_se->second[e].second,eltNum);
2354 }
2355 }
2356 }
2357
2358 //! construct a 2d/3d extruded mesh from a 1d/2d section mesh with quadrangle in 2D, with prism in 3D
2359 /*
2360 ms : mesh of the section
2361 O : section origin
2362 D : section end (direction of extrusion is OD)
2363 nbl : number of layers
2364 namingDomain, namingSection, namingSide : 0 means no naming, 1 means unique name, >1 means each intern domain/section/side is named
2365
2366 domain names are inherited from section domain names
2367
2368 for instance, assume the following mesh of a rectangle
2369 Gamma
2370 -----------------
2371 | | |
2372 SigmaP | Omg1 | Omg2 | SigmaM
2373 | | |
2374 -----------------
2375 Gamma
2376
2377 namingDomain=0 will produce only the full 3d domain "Omega"
2378 namingDomain=1 will produce the volumic domains : "Omg1_e", "Omg2_e"
2379 namingDomain=2 will produce the volumic domains : "Omg1_e1","Omg1_e2", ..., "Omg1_en", "Omg2_e1","Omg2_e2",..., "Omg2_en"
2380 namingSection=0 will not produce any section domains
2381 namingSection=1 will produce the section domains : "Omg1_0", "Omg2_0" (section 0), "Omg1_n", "Omg2_n" (section n)
2382 namingSection=2 will produce the section domains : "Omg1_0", "Omg2_0" (section 0), "Omg1_1", "Omg2_2" (section 2),..., "Omg1_n", "Omg2_n" (section n)
2383 namingSide=0 will not produce any side domains
2384 namingSide=1 will produce the side domains : "Gamma_e", "SigmaP_e, "SigmaM_e"
2385 namingSide=2 will produce the side domains : "Gamma_e1","Gamma_e2", ..., "Gamma_en", "SigmaM_e1", "SigmaM_e2", ..., "SigmaM_en", "SigmaP_e1", "SigmaP_e2", ..., "SigmaP_en"
2386
2387 Default name for volumic domain is "Omega", for section domain is "Section"
2388 If section has no side domain defined, the extruded mesh will not have any side domains
2389 SideOfSide domain in 3D are not yet managed
2390
2391 Remarks :
2392 Only one order mesh are currently extruded
2393 nbelt_section * n extruded elements are ordered by layer
2394
2395
2396 */
buildExtrusion(const Mesh & ms,const Point & O,const Point & D,number_t nbl,number_t namingDomain,number_t namingSection,number_t namingSide)2397 void Mesh::buildExtrusion(const Mesh& ms, const Point& O, const Point& D, number_t nbl,
2398 number_t namingDomain, number_t namingSection, number_t namingSide)
2399 {
2400 trace_p->push("Mesh::(Mesh,...) extrusion");
2401 dimen_t dim_mesh=ms.meshDim();
2402 if(dim_mesh!=1 && dim_mesh!=2) error("dim_not_in_range", 1, 2);
2403
2404 if(name_=="") name_=ms.name_+"_extrusion";
2405 comment_="extrusion of "+ms.name_+" using "+tostring(nbl)+" layers";
2406 isMadeOfSimplices_=false;
2407 order_=1;
2408 lastIndex_=0;
2409 firstOrderMesh_p=this;
2410 geometry_p=0;
2411 number_t nbelt_s=ms.nbOfElements();
2412 number_t nbver_s=ms.nbOfVertices();
2413 elements_.resize(nbelt_s*nbl);
2414 nodes.resize(nbver_s*(nbl+1));
2415 dimen_t dim_space=ms.spaceDim(); //space dimension of nodes of the section mesh
2416 Point dir=D-O;
2417 real_t dt=norm2(dir)/nbl;
2418
2419 //create nodes (only from vertices), ordered by section
2420 // M_k = M + O + k*dt*OD k=0, nbl with ds= ||OD||/nbl
2421 dimen_t dim=std::max(std::max(dimen_t(dim_mesh+1),dim_space),dir.dim()); //final dimension of nodes
2422 Point Oe(O,dim), dire(dir,dim);
2423 number_t m=0;
2424 std::vector<number_t>::const_iterator itv;
2425 for(number_t k=0; k<=nbl; k++)
2426 {
2427 Point Q=Oe+k*dt*dire;
2428 for(itv=ms.vertices_.begin(); itv!=ms.vertices_.end(); ++itv, m++)
2429 nodes[m]=Point(ms.nodes[*itv-1],dim)+Q;
2430 }
2431
2432 //construct vertex indices (begins at 1)
2433 vertices_.resize(nodes.size());
2434 for(number_t i = 0; i < nodes.size(); i++) { vertices_[i] = i + 1; }
2435
2436
2437 //create elements (ordered by section)
2438 ShapeType sh=_segment;
2439 Interpolation* interp_p = findInterpolation(Lagrange, standard, 1, H1);
2440 RefElement* ref_p = findRefElement(_quadrangle, interp_p);
2441 if(dim_mesh==2) //2D case, assuming same element shapetype
2442 {
2443 sh=(*ms.elements_.begin())->shapeType();
2444 if(sh==_triangle) ref_p= findRefElement(_prism, interp_p);
2445 else if(sh==_quadrangle) ref_p= findRefElement(_hexahedron, interp_p);
2446 else error("shape_not_handled", words("shape", sh));
2447 }
2448 number_t e=0;
2449 for(number_t k=0; k<nbl; k++) //loop on slices
2450 {
2451 std::vector<GeomElement*>::const_iterator ites=ms.elements_.begin();
2452 for(; ites!=ms.elements_.end(); ++ites, e++)
2453 {
2454 elements_[e]=new GeomElement(this, ref_p, dim, e+1);
2455 MeshElement* melt=elements_[e]->meshElement();
2456 number_t nv=(*ites)->numberOfVertices();
2457 melt->vertexNumbers.resize(2*nv);
2458 if(sh==_segment) //1D case
2459 {
2460 melt->vertexNumbers[0]=(*ites)->vertexNumber(1)+k*nbver_s;
2461 melt->vertexNumbers[1]=(*ites)->vertexNumber(2)+k*nbver_s;
2462 melt->vertexNumbers[2]=(*ites)->vertexNumber(2)+(k+1)*nbver_s;
2463 melt->vertexNumbers[3]=(*ites)->vertexNumber(1)+(k+1)*nbver_s;
2464 }
2465 else
2466 {
2467 for(number_t i=0; i<nv; ++i)
2468 {
2469 melt->vertexNumbers[i]=(*ites)->vertexNumber(i+1)+k*nbver_s;
2470 melt->vertexNumbers[nv+i]=melt->vertexNumbers[i]+nbver_s;
2471 }
2472 }
2473 melt->nodeNumbers = melt->vertexNumbers;
2474 melt->setNodes(nodes);
2475 }
2476 }
2477
2478 // construct vertex indices (begins at 1)
2479 vertices_.resize(nodes.size());
2480 for(number_t i = 0; i < nodes.size(); i++) { vertices_[i] = i + 1; }
2481
2482 //create domains
2483 std::vector<GeomDomain*>::const_iterator itd=ms.domains().begin();
2484 std::vector<GeomDomain*> doms;
2485 for(; itd!=ms.domains().end(); ++itd)
2486 if((*itd)->dim()==dim_mesh && (*itd)->domType()==_meshDomain) doms.push_back(*itd);
2487
2488 if(namingDomain <2) //unique domain
2489 {
2490 if(doms.size()==0 || namingDomain==0) //only the full domain named Omega
2491 {
2492 MeshDomain* meshdom_p = (new GeomDomain(*this, "Omega", dim_mesh+1))->meshDomain();
2493 meshdom_p->geomElements = elements_;
2494 domains_.push_back(meshdom_p);
2495 }
2496 else
2497 for(number_t d=0; d<doms.size(); d++)
2498 {
2499 number_t e=0;
2500 MeshDomain* meshdom_p = (new GeomDomain(*this, doms[d]->name()+"_e", dim_mesh+1))->meshDomain();
2501 domains_.push_back(meshdom_p);
2502 MeshDomain* dom=doms[d]->meshDomain();
2503 number_t nbe = dom->numberOfElements();
2504 meshdom_p->geomElements.resize(nbe*nbl);
2505 for(number_t k=0; k<nbl; k++)
2506 for(number_t i=0; i<nbe; i++, e++)
2507 meshdom_p->geomElements[e]=elements_[k*nbelt_s+dom->geomElements[i]->number()-1];
2508 }
2509 }
2510 if(namingDomain >= 2) //one domain by slice
2511 {
2512 if(doms.size()==0)
2513 {
2514 number_t e=0;
2515 for(number_t k=0; k<nbl; k++) //loop on slices
2516 {
2517 MeshDomain* meshdom_p = (new GeomDomain(*this, "Omega_e"+tostring(k+1), dim_mesh+1))->meshDomain();
2518 meshdom_p->geomElements.resize(nbelt_s);
2519 for(number_t i=0; i<nbelt_s; i++, e++) meshdom_p->geomElements[i]=elements_[e];
2520 domains_.push_back(meshdom_p);
2521 }
2522 }
2523 else
2524 for(number_t d=0; d<doms.size(); d++)
2525 {
2526 MeshDomain* dom=doms[d]->meshDomain();
2527 number_t nbe = dom->numberOfElements();
2528 string_t dnam = doms[d]->name();
2529 for(number_t k=0; k<nbl; k++)
2530 {
2531 number_t e=0;
2532 MeshDomain* meshdom_p = (new GeomDomain(*this, dnam+"_e"+tostring(k+1), dim_mesh+1))->meshDomain();
2533 domains_.push_back(meshdom_p);
2534 meshdom_p->geomElements.resize(nbe);
2535 for(number_t i=0; i<nbe; i++, e++)
2536 meshdom_p->geomElements[e]=elements_[k*nbe+dom->geomElements[i]->number()-1];
2537 }
2538 }
2539 }
2540
2541 number_t last = elements_.size()+1;
2542 if(namingSection >0) //create section domain
2543 {
2544 number_t s1=1;
2545 if(sh==_quadrangle) s1=6; //left side numbering
2546 number_t s2=3;
2547 if(sh==_triangle) s2=5; //right side numbering
2548 if(namingSection==1)
2549 {
2550 if(doms.size()==0)
2551 {
2552 MeshDomain* meshdom_p = (new GeomDomain(*this, "Section_0", dim_mesh))->meshDomain();
2553 domains_.push_back(meshdom_p);
2554 meshdom_p->geomElements.resize(nbelt_s);
2555 for(number_t i=0; i<nbelt_s; i++, last++)
2556 meshdom_p->geomElements[i]=new GeomElement(elements_[i],s1,last);
2557 meshdom_p = (new GeomDomain(*this, "Section_"+tostring(nbl), dim_mesh))->meshDomain();
2558 domains_.push_back(meshdom_p);
2559 meshdom_p->geomElements.resize(nbelt_s);
2560 for(number_t i=0; i<nbelt_s; i++, last++)
2561 meshdom_p->geomElements[i]=new GeomElement(elements_[(nbl-1)*nbelt_s+i],s2,last);
2562 }
2563 else
2564 {
2565 for(number_t d=0; d<doms.size(); d++)
2566 {
2567 MeshDomain* dom=doms[d]->meshDomain();
2568 number_t nbe = dom->numberOfElements();
2569 string_t dnam = doms[d]->name();
2570 MeshDomain* meshdom_p = (new GeomDomain(*this, dnam+"_0", dim_mesh))->meshDomain();
2571 domains_.push_back(meshdom_p);
2572 meshdom_p->geomElements.resize(nbe);
2573 for(number_t i=0; i<nbe; i++, last++)
2574 meshdom_p->geomElements[i]=new GeomElement(elements_[dom->geomElements[i]->number()-1],s1,last);
2575 meshdom_p = (new GeomDomain(*this, dnam+"_"+tostring(nbl), dim_mesh))->meshDomain();
2576 domains_.push_back(meshdom_p);
2577 meshdom_p->geomElements.resize(nbe);
2578 for(number_t i=0; i<nbe; i++, last++)
2579 meshdom_p->geomElements[i]=new GeomElement(elements_[(nbl-1)*nbelt_s+dom->geomElements[i]->number()-1],s2,last);
2580 }
2581 }
2582 }
2583 else
2584 {
2585 if(doms.size()==0)
2586 {
2587 for(number_t k=0; k<nbl+1; k++)
2588 {
2589 MeshDomain* meshdom_p = (new GeomDomain(*this, "Section_"+tostring(k), dim_mesh))->meshDomain();
2590 domains_.push_back(meshdom_p);
2591 meshdom_p->geomElements.resize(nbelt_s);
2592 for(number_t i=0; i<nbelt_s; i++, last++)
2593 {
2594 if(k<nbl) meshdom_p->geomElements[i]=new GeomElement(elements_[k*nbelt_s+i],s1,last);
2595 else meshdom_p->geomElements[i]=new GeomElement(elements_[(k-1)*nbelt_s+i],s2,last);
2596 if(k>0 && k<nbl) meshdom_p->geomElements[i]->parentSides().push_back(GeoNumPair(elements_[(k-1)*nbelt_s+i], s2));
2597 }
2598 }
2599 }
2600 else
2601 {
2602 for(number_t d=0; d<doms.size(); d++)
2603 {
2604 MeshDomain* dom=doms[d]->meshDomain();
2605 number_t nbe = dom->numberOfElements();
2606 string_t dnam = doms[d]->name();
2607 for(number_t k=0; k<nbl+1; k++)
2608 {
2609 MeshDomain* meshdom_p = (new GeomDomain(*this, dnam+"_"+tostring(k), dim_mesh))->meshDomain();
2610 domains_.push_back(meshdom_p);
2611 meshdom_p->geomElements.resize(nbe);
2612 for(number_t i=0; i<nbe; i++, last++)
2613 {
2614 number_t j= dom->geomElements[i]->number()-1;
2615 if(k<nbl) meshdom_p->geomElements[i]=new GeomElement(elements_[k*nbelt_s+j],s1,last);
2616 else meshdom_p->geomElements[i]=new GeomElement(elements_[(k-1)*nbelt_s+j],s2,last);
2617 if(k>0 && k<nbl) meshdom_p->geomElements[i]->parentSides().push_back(GeoNumPair(elements_[(k-1)*nbelt_s+j], s2));
2618 }
2619 }
2620 }
2621 }
2622 }
2623 }
2624
2625 if(namingSide >0) //create side domain, propagates sidename of section mesh
2626 {
2627 itd=ms.domains().begin();
2628 doms.clear();
2629 for(; itd!=ms.domains().end(); ++itd)
2630 if((*itd)->dim()==dim_mesh-1 && (*itd)->domType()==_meshDomain) doms.push_back(*itd);
2631
2632 std::vector<number_t> sToS; // side numbering of section element to side numbering of side of extruded element
2633 switch(sh)
2634 {
2635 case _segment : sToS.resize(3,0); sToS[1]=4; sToS[2]=2; break;
2636 case _triangle : sToS.resize(4,0); sToS[1]=2; sToS[2]=3; sToS[3]=4; break;
2637 case _quadrangle: sToS.resize(5,0); sToS[1]=1; sToS[2]=2; sToS[3]=4; sToS[4]=5; break;
2638 default : error("shape_not_handled", words("shape", sh));
2639 }
2640 if(doms.size()>0)
2641 {
2642 for(itd=doms.begin(); itd!=doms.end(); ++itd)
2643 {
2644 MeshDomain* sidom=(*itd)->meshDomain();
2645 number_t nbe = sidom->numberOfElements();
2646 if(namingSide==1) //only one boundary
2647 {
2648 MeshDomain* meshdom_p = (new GeomDomain(*this, (*itd)->name()+"_e", dim_mesh))->meshDomain();
2649 domains_.push_back(meshdom_p);
2650 meshdom_p->geomElements.resize(nbe*nbl);
2651 number_t e=0;
2652 for(number_t k=0; k<nbl; k++) //loop on slice
2653 {
2654 for(number_t i=0; i<nbe; i++, last++, e++)
2655 {
2656 number_t ind = sidom->geomElements[i]->parentSides()[0].first->number();
2657 number_t sid = sToS[sidom->geomElements[i]->parentSides()[0].second];
2658 meshdom_p->geomElements[e]=new GeomElement(elements_[k*nbelt_s+ind-1],sid,last);
2659 }
2660 }
2661 }
2662 else //one boundary by layer
2663 {
2664 for(number_t k=0; k<nbl; k++) //loop on slice
2665 {
2666 MeshDomain* meshdom_p = (new GeomDomain(*this, (*itd)->name()+"_e"+tostring(k), dim_mesh))->meshDomain();
2667 domains_.push_back(meshdom_p);
2668 meshdom_p->geomElements.resize(nbe);
2669 for(number_t i=0; i<nbe; i++, last++)
2670 {
2671 number_t ind = sidom->geomElements[i]->parentSides()[0].first->number();
2672 number_t sid = sToS[sidom->geomElements[i]->parentSides()[0].second];
2673 meshdom_p->geomElements[i]=new GeomElement(elements_[k*nbelt_s+ind-1],sid,last);
2674 }
2675 }
2676 }
2677 }
2678 }
2679 }
2680
2681 // compute measures and orientation of mesh elements
2682 buildGeomData();
2683 setShapeTypes();
2684 trace_p->pop();
2685 }
2686
2687 /*create pyramid mesh from hexaedron mesh
2688 by spliting each hexaedron into six pyramids using the hexahedron centroids as new nodes c
2689 hexahedron face 1: 2 6 5 1 -> pyramid 2 6 5 1 c
2690 hexahedron face 2: 7 6 2 3 -> pyramid 7 6 2 3 c
2691 hexahedron face 3: 5 6 7 8 -> pyramid 5 6 7 8 c
2692 hexahedron face 4: 3 4 8 7 -> pyramid 3 4 8 7 c
2693 hexahedron face 5: 8 4 1 5 -> pyramid 8 4 1 5 c
2694 hexahedron face 6: 1 4 3 2 -> pyramid 1 4 3 2 c
2695 split only first order hexaedra
2696 assume current void mesh
2697 */
buildPyramidFromHexadron(const Mesh & hexMesh)2698 void Mesh::buildPyramidFromHexadron(const Mesh& hexMesh)
2699 {
2700 trace_p->push("Mesh::buildPyramidFromHexadron");
2701
2702 if(name_=="") name_=hexMesh.name_+"_toPyramids";
2703 comment_="mesh from spliting hexadron into six pyramids";
2704 isMadeOfSimplices_=false;
2705 order_=1;
2706 lastIndex_=0;
2707 firstOrderMesh_p=this;
2708 geometry_p=0;
2709 number_t nbelt_h=hexMesh.nbOfElements();
2710 number_t nbver_h=hexMesh.nbOfVertices();
2711 elements_.resize(6*nbelt_h);
2712 nodes.resize(nbver_h+nbelt_h);
2713
2714 //create nodes only from vertices
2715 std::vector<number_t>::const_iterator itv;
2716 std::vector<Point>::iterator itn=nodes.begin();
2717 for(itv=hexMesh.vertices_.begin(); itv!=hexMesh.vertices_.end(); ++itv, ++itn) //copy original nodes
2718 *itn=hexMesh.nodes[*itv-1];
2719 std::vector<GeomElement*>::const_iterator iteh=hexMesh.elements_.begin();
2720 for(; iteh!=hexMesh.elements_.end(); ++iteh, ++itn)
2721 {
2722 if((*iteh)->shapeType()!=_hexahedron)
2723 {
2724 where("Mesh::buildPyramidFromHexadron()");
2725 error("shape_not_handled", words("shape", (*iteh)->shapeType()));
2726 }
2727 *itn = (*iteh)->meshElement()->centroid;
2728 }
2729
2730 //construct vertex indices (begins at 1)
2731 vertices_.resize(nodes.size());
2732 for(number_t i = 0; i < nodes.size(); i++) { vertices_[i] = i + 1; }
2733
2734 //create elements
2735 std::vector<GeomElement*>::iterator ite=elements_.begin();
2736 number_t nc = nbver_h+1, e=0;
2737 Interpolation* interp_p = findInterpolation(Lagrange, standard, 1, H1);
2738 RefElement* ref_p = findRefElement(_pyramid, interp_p);
2739 const GeomRefElement* gref_h = findRefElement(_hexahedron, interp_p)->geomRefElement();
2740 for(iteh=hexMesh.elements_.begin(); iteh!=hexMesh.elements_.end(); ++iteh, nc++)
2741 {
2742 for(number_t s=1; s<=6; s++, e++, ++ite) //loop on hexahedron faces
2743 {
2744 *ite = new GeomElement(this, ref_p, 3, e+1);
2745 MeshElement* melt=(*ite)->meshElement();
2746 melt->vertexNumbers.resize(5);
2747 for(number_t i=1; i<=4; ++i)
2748 melt->vertexNumbers[i-1]=(*iteh)->vertexNumber(gref_h->sideVertexNumber(i,s));
2749 melt->vertexNumbers[4] = nc;
2750 melt->nodeNumbers = melt->vertexNumbers;
2751 melt->setNodes(nodes);
2752 }
2753 }
2754 number_t last =elements_.size()+1;
2755
2756 //create domains with same name
2757 std::vector<GeomDomain*>::const_iterator itdh=hexMesh.domains().begin();
2758 for(; itdh!=hexMesh.domains().end(); ++itdh)
2759 {
2760 MeshDomain* domh=(*itdh)->meshDomain();
2761 number_t nbe = domh->numberOfElements();
2762 if((*itdh)->dim()==3) //volumic domain
2763 {
2764 MeshDomain* domp = (new GeomDomain(*this, (*itdh)->name(), 3))->meshDomain();
2765 domains_.push_back(domp);
2766 domp->geomElements.resize(6*nbe);
2767 ite=domp->geomElements.begin();
2768 for(iteh = domh->geomElements.begin(); iteh!=domh->geomElements.end(); ++iteh)
2769 {
2770 number_t num=(*iteh)->number()-1;
2771 for(number_t k=0; k<6; ++k,++ite) *ite = elements_[6*num+k];
2772 }
2773 }
2774 if((*itdh)->dim()==2) //surfacic domain
2775 {
2776 MeshDomain* domp = (new GeomDomain(*this, (*itdh)->name(), 2))->meshDomain();
2777 domains_.push_back(domp);
2778 domp->geomElements.resize(nbe);
2779 ite=domp->geomElements.begin();
2780 for(iteh = domh->geomElements.begin(); iteh!=domh->geomElements.end(); ++iteh, ++ite)
2781 {
2782 const std::vector<GeoNumPair>& psh= (*iteh)->parentSides();
2783 std::vector<GeoNumPair>::const_iterator itsh=psh.begin();
2784 GeomElement* parp = elements_[6*(itsh->first->number()-1)+itsh->second-1];
2785 GeomElement* geop = new GeomElement(parp,1,last++);
2786 *ite = geop;
2787 itsh++;
2788 for(; itsh!=psh.end(); ++itsh)
2789 {
2790 parp =elements_[6*(itsh->first->number()-1)+itsh->second-1];
2791 geop->parentSides().push_back(GeoNumPair(parp,1));
2792 }
2793 }
2794 }
2795 }
2796
2797 // compute measures and orientation of mesh elements
2798 buildGeomData();
2799 setShapeTypes();
2800 trace_p->pop();
2801 }
2802
2803 } // end of namespace xlifepp
2804