v0.14.0
Loading...
Searching...
No Matches
RefElementMultiIndices.cpp
Go to the documentation of this file.
1/** \file RefElementMultiIndices.cpp
2 * \brief Multi-index containers for finite elements
3 */
4
5
6namespace MoFEM {
7
8// ref moab ent
10
11const boost::shared_ptr<SideNumber> RefElement::nullSideNumber =
12 boost::shared_ptr<SideNumber>();
13
14// ref moab FiniteElement
15RefElement::RefElement(const boost::shared_ptr<RefEntity> &ref_ents_ptr)
16 : interface_RefEntity<RefEntity>(ref_ents_ptr) {}
17
18std::ostream &operator<<(std::ostream &os, const RefElement &e) {
19 os << " ref egdes " << e.getBitRefEdges();
20 os << " " << *(e.sPtr);
21 return os;
22}
23
25 const boost::shared_ptr<RefEntity> &ref_ents_ptr)
26 : RefElement(ref_ents_ptr) {
27 switch (ref_ents_ptr->getEntType()) {
28 case MBENTITYSET:
29 break;
30 default:
31 THROW_MESSAGE("this work only for MESHSETs");
32 }
33}
34const boost::shared_ptr<SideNumber> &
36 NOT_USED(ent);
37 SideNumber_multiIndex::iterator miit;
38 miit =
40 .insert(boost::shared_ptr<SideNumber>(new SideNumber(ent, 0, 0, 0)))
41 .first;
42 return *miit;
43}
45 const boost::shared_ptr<RefEntity> &ref_ents_ptr)
46 : RefElement(ref_ents_ptr) {
47 Tag th_RefBitEdge;
48 moab::Interface &moab = getRefEntityPtr()->getBasicDataPtr()->moab;
49 rval = moab.tag_get_handle("_RefBitEdge", th_RefBitEdge);
51 rval = moab.tag_get_by_ptr(th_RefBitEdge, &ref_ents_ptr->ent, 1,
52 (const void **)&tag_BitRefEdges);
54 switch (ref_ents_ptr->getEntType()) {
55 case MBPRISM:
56 break;
57 default:
58 THROW_MESSAGE("this work only for PRISMs");
59 }
60 EntityHandle prism = getEnt();
61 int num_nodes;
62 const EntityHandle *conn;
63 rval = moab.get_connectivity(prism, conn, num_nodes, true);
65 assert(num_nodes == 6);
66 for (int nn = 0; nn != 6; ++nn) {
68 .insert(
69 boost::shared_ptr<SideNumber>(new SideNumber(conn[nn], nn, 0, -1)));
70 }
71 // Range face_side3, face_side4;
72 // CHKERR moab.get_adjacencies(conn, 3, 2, true, face_side3);
73 // CHKERR moab.get_adjacencies(&conn[3], 3, 2, true, face_side4);
74 // if (face_side3.size() != 1)
75 // THROW_MESSAGE("prism don't have side face 3");
76 // if (face_side4.size() != 1)
77 // THROW_MESSAGE("prims don't have side face 4");
78 // getSideNumberPtr(*face_side3.begin());
79 // getSideNumberPtr(*face_side4.begin());
80}
81const boost::shared_ptr<SideNumber> &
83
84 moab::Interface &moab = getRefEntityPtr()->getBasicDataPtr()->moab;
85
86 SideNumber_multiIndex::iterator miit = side_number_table.find(ent);
87 // this int is in table then return pointer
88 if (miit != side_number_table.end())
89 return *miit;
90
91 // if ent is a this prism
92 if (sPtr->ent == ent) {
93 miit =
95 .insert(boost::shared_ptr<SideNumber>(new SideNumber(ent, 0, 0, 0)))
96 .first;
97 return *miit;
98 }
99
100 // if ent is meshset
101 if (type_from_handle(ent) == MBENTITYSET) {
102 miit =
104 .insert(boost::shared_ptr<SideNumber>(new SideNumber(ent, 0, 0, 0)))
105 .first;
106 return *miit;
107 }
108
109 // use moab to get sense, side and offset
110
111 int side_number, sense, offset;
112 rval = moab.side_number(sPtr->ent, ent, side_number, sense, offset);
113
114 // it has to be degenerated prism, get sense from nodes topology
115 if (side_number == -1 || rval != MB_SUCCESS) {
116
117 if (type_from_handle(ent) == MBVERTEX) {
118 THROW_MESSAGE("Huston we have problem, vertex (specified by ent) is not "
119 "part of prism, that is impossible (top tip: check your "
120 "prisms)");
121 }
122
123 // get prism connectivity
124 int num_nodes;
125 const EntityHandle *conn;
126 rval = moab.get_connectivity(sPtr->ent, conn, num_nodes, true);
128 assert(num_nodes == 6);
129 // get ent connectivity
130 const EntityHandle *conn_ent;
131 rval = moab.get_connectivity(ent, conn_ent, num_nodes, true);
133
134 // for(int nn = 0; nn<6;nn++) {
135 // std::cerr << conn[nn] << " ";
136 // };
137 // std::cerr << std::endl;
138 // for(int nn = 0; nn<num_nodes;nn++) {
139 // std::cerr << conn_ent[nn] << " ";
140 // }
141 // std::cerr << std::endl;
142
143 // bottom face
144 EntityHandle face3[3] = {conn[0], conn[1], conn[2]};
145 // top face
146 EntityHandle face4[3] = {conn[3], conn[4], conn[5]};
147 if (num_nodes == 3) {
148 int sense_p1_map[3][3] = {{0, 1, 2}, {1, 2, 0}, {2, 0, 1}};
149 int sense_m1_map[3][3] = {{0, 2, 1}, {1, 0, 2}, {2, 1, 0}};
150 EntityHandle *conn0_3_ptr = std::find(face3, &face3[3], conn_ent[0]);
151 if (conn0_3_ptr != &face3[3]) {
152 offset = std::distance(face3, conn0_3_ptr);
153 if (face3[sense_p1_map[offset][0]] == conn_ent[0] &&
154 face3[sense_p1_map[offset][1]] == conn_ent[1] &&
155 face3[sense_p1_map[offset][2]] == conn_ent[2]) {
156 miit = const_cast<SideNumber_multiIndex &>(side_number_table)
157 .insert(boost::shared_ptr<SideNumber>(
158 new SideNumber(ent, 3, 1, offset)))
159 .first;
160 return *miit;
161 } else if (face3[sense_m1_map[offset][0]] == conn_ent[0] &&
162 face3[sense_m1_map[offset][1]] == conn_ent[1] &&
163 face3[sense_m1_map[offset][2]] == conn_ent[2]) {
164 miit = const_cast<SideNumber_multiIndex &>(side_number_table)
165 .insert(boost::shared_ptr<SideNumber>(
166 new SideNumber(ent, 3, -1, offset)))
167 .first;
168 return *miit;
169 }
170 }
171 EntityHandle *conn0_4_ptr = std::find(face4, &face4[3], conn_ent[0]);
172 if (conn0_4_ptr != &face4[3]) {
173 offset = std::distance(face4, conn0_4_ptr);
174 if (face4[sense_p1_map[offset][0]] == conn_ent[0] &&
175 face4[sense_p1_map[offset][1]] == conn_ent[1] &&
176 face4[sense_p1_map[offset][2]] == conn_ent[2]) {
177 miit = const_cast<SideNumber_multiIndex &>(side_number_table)
178 .insert(boost::shared_ptr<SideNumber>(
179 new SideNumber(ent, 4, 1, 3 + offset)))
180 .first;
181 return *miit;
182 } else if (face4[sense_m1_map[offset][0]] == conn_ent[0] &&
183 face4[sense_m1_map[offset][1]] == conn_ent[1] &&
184 face4[sense_m1_map[offset][2]] == conn_ent[2]) {
185 miit = const_cast<SideNumber_multiIndex &>(side_number_table)
186 .insert(boost::shared_ptr<SideNumber>(
187 new SideNumber(ent, 4, -1, 3 + offset)))
188 .first;
189 return *miit;
190 } else {
191 std::cerr << conn_ent[0] << " " << conn_ent[1] << " " << conn_ent[2]
192 << std::endl;
193 std::cerr << face3[0] << " " << face3[1] << " " << face3[2]
194 << std::endl;
195 std::cerr << face4[0] << " " << face4[1] << " " << face4[2]
196 << std::endl;
197 std::cerr << offset << std::endl;
198 THROW_MESSAGE("Huston we have problem");
199 }
200 }
201 THROW_MESSAGE("Huston we have problem");
202 }
203
204 if (num_nodes == 2) {
205 {
206 // Triangle edges
207 EntityHandle edges[6][2] = {
208 {conn[0], conn[1]} /*0*/, {conn[1], conn[2]} /*1*/,
209 {conn[2], conn[0]} /*2*/, {conn[3], conn[4]} /*3+3*/,
210 {conn[4], conn[5]} /*3+4*/, {conn[5], conn[3]} /*3+5*/
211 };
212 for (int ee = 0; ee < 6; ee++) {
213 if (((conn_ent[0] == edges[ee][0]) &&
214 (conn_ent[1] == edges[ee][1])) ||
215 ((conn_ent[0] == edges[ee][1]) &&
216 (conn_ent[1] == edges[ee][0]))) {
217 side_number = ee;
218 if (ee >= 3) {
219 side_number += 3;
220 EntityHandle *conn0_4_ptr =
221 std::find(face4, &face4[3], conn_ent[0]);
222 offset = std::distance(face4, conn0_4_ptr) + 3;
223 } else {
224 EntityHandle *conn0_3_ptr =
225 std::find(face3, &face3[3], conn_ent[0]);
226 offset = std::distance(face3, conn0_3_ptr);
227 }
228 sense = 1;
229 if ((conn_ent[0] == edges[ee][1]) && (conn_ent[1] == edges[ee][0]))
230 sense = -1;
231 miit = const_cast<SideNumber_multiIndex &>(side_number_table)
232 .insert(boost::shared_ptr<SideNumber>(
233 new SideNumber(ent, side_number, sense, offset)))
234 .first;
235 return *miit;
236 }
237 }
238 }
239 // {
240 // // Edges through thickness
241 // EntityHandle edges[3][2] = {
242 // { conn[0], conn[3] }, { conn[1], conn[4] }, { conn[2], conn[5] }
243 // };
244 // for(int ee = 0;ee<3;ee++) {
245 // if(
246 // (( conn_ent[0] == edges[ee][0] )&&( conn_ent[1] == edges[ee][1]
247 // ))||
248 // (( conn_ent[0] == edges[ee][1] )&&( conn_ent[1] == edges[ee][0]
249 // ))
250 // ) {
251 // side_number = 3+ee;
252 // offset = std::distance(conn,find(conn,&conn[6],conn_ent[0]));
253 // sense = 1;
254 // if(( conn_ent[0] == edges[ee][1] )&&( conn_ent[1] == edges[ee][0]
255 // )) sense = -1; miit =
256 // const_cast<SideNumber_multiIndex&>(side_number_table).insert(SideNumber(ent,side_number,sense,offset)).first;
257 // return const_cast<SideNumber*>(&*miit);
258 // }
259 // }
260 // }
261 // for(int nn = 0; nn<6;nn++) {
262 // std::cerr << conn[nn] << " ";
263 // };
264 // std::cerr << std::endl;
265 // std::cerr << conn_ent[0] << " " << conn_ent[1] << std::endl;
266 THROW_MESSAGE("Huston we have problem");
267 }
268 std::ostringstream sss;
269 sss << "this not working: " << ent << " type: " << type_from_handle(ent)
270 << " " << MBEDGE << " " << MBTRI << std::endl;
271 THROW_MESSAGE(sss.str().c_str());
272 }
273 miit = const_cast<SideNumber_multiIndex &>(side_number_table)
274 .insert(boost::shared_ptr<SideNumber>(
275 new SideNumber(ent, side_number, sense, offset)))
276 .first;
277 return *miit;
278 THROW_MESSAGE("not implemented");
279 return nullSideNumber;
280}
281
282RefElementVolume::RefElementVolume(const boost::shared_ptr<RefEntity> &ref_ents_ptr)
283 : RefElement(ref_ents_ptr), tag_BitRefEdges(NULL) {
284 Tag th_RefBitEdge;
285 moab::Interface &moab = getRefEntityPtr()->getBasicDataPtr()->moab;
286 rval = moab.tag_get_handle("_RefBitEdge", th_RefBitEdge);
288 rval = moab.tag_get_by_ptr(th_RefBitEdge, &ref_ents_ptr->ent, 1,
289 (const void **)&tag_BitRefEdges);
291 switch (ref_ents_ptr->getEntType()) {
292 case MBTET:
293 case MBHEX:
294 break;
295 default:
296 THROW_MESSAGE("this work only for TETs or HEXs");
297 }
299 .insert(boost::make_shared<SideNumber>(sPtr->ent, 0, 0, 0));
300}
301
302const boost::shared_ptr<SideNumber> &
304 moab::Interface &moab = getRefEntityPtr()->getBasicDataPtr()->moab;
305 auto miit = side_number_table.find(ent);
306 if (miit != side_number_table.end())
307 return *miit;
308 if (sPtr->ent == ent) {
309 miit = const_cast<SideNumber_multiIndex &>(side_number_table)
310 .insert(boost::make_shared<SideNumber>(ent, 0, 0, 0))
311 .first;
312 return *miit;
313 }
314 if (type_from_handle(ent) == MBENTITYSET) {
315 miit = const_cast<SideNumber_multiIndex &>(side_number_table)
316 .insert(boost::make_shared<SideNumber>(ent, 0, 0, 0))
317 .first;
318 return *miit;
319 }
320 int side_number, sense, offset;
321 rval = moab.side_number(sPtr->ent, ent, side_number, sense, offset);
323 auto p_miit = const_cast<SideNumber_multiIndex &>(side_number_table)
324 .insert(boost::make_shared<SideNumber>(ent, side_number,
325 sense, offset));
326 miit = p_miit.first;
327 if (miit->get()->ent != ent)
328 THROW_MESSAGE("this not working");
329
330 return *miit;
331}
332std::ostream &operator<<(std::ostream &os, const RefElementVolume &e) {
333 os << " ref egdes " << e.getBitRefEdges();
334 os << " " << *e.sPtr;
335 return os;
336}
337
338RefElementFace::RefElementFace(const boost::shared_ptr<RefEntity> &ref_ents_ptr)
339 : RefElement(ref_ents_ptr) {
340
341 int nb_nodes = 0;
342 int nb_edges = 0;
343 switch (ref_ents_ptr->getEntType()) {
344 case MBTRI:
345 nb_nodes = nb_edges = 3;
346 break;
347 case MBQUAD:
348 nb_nodes = nb_edges = 4;
349 break;
350 default:
351 THROW_MESSAGE("this works only for TRIs and QUADs");
352 }
353 int side_number, sense, offset;
354 EntityHandle tri = getEnt();
355 int num_nodes;
356 const EntityHandle *conn;
357 moab::Interface &moab = getRefEntityPtr()->getBasicDataPtr()->moab;
358 rval = moab.get_connectivity(tri, conn, num_nodes, true);
360 for (int nn = 0; nn < nb_nodes; nn++) {
362 .insert(
363 boost::shared_ptr<SideNumber>(new SideNumber(conn[nn], nn, 0, 0)));
364 }
365 for (int ee = 0; ee < nb_edges; ee++) {
366 EntityHandle edge;
367 rval = moab.side_element(tri, 1, ee, edge);
369 rval = moab.side_number(tri, edge, side_number, sense, offset);
372 .insert(boost::shared_ptr<SideNumber>(
373 new SideNumber(edge, ee, sense, offset)));
374 }
376 .insert(boost::shared_ptr<SideNumber>(new SideNumber(tri, 0, 0, 0)));
377}
378const boost::shared_ptr<SideNumber> &
380 moab::Interface &moab = getRefEntityPtr()->getBasicDataPtr()->moab;
381 SideNumber_multiIndex::iterator miit = side_number_table.find(ent);
382 if (miit != side_number_table.end())
383 return *miit;
384 if (sPtr->ent == ent) {
385 miit =
387 .insert(boost::shared_ptr<SideNumber>(new SideNumber(ent, 0, 0, 0)))
388 .first;
389 return *miit;
390 }
391 if (type_from_handle(ent) == MBENTITYSET) {
392 miit =
394 .insert(boost::shared_ptr<SideNumber>(new SideNumber(ent, 0, 0, 0)))
395 .first;
396 return *miit;
397 }
398 int side_number, sense, offset;
399 rval = moab.side_number(sPtr->ent, ent, side_number, sense, offset);
401 miit = const_cast<SideNumber_multiIndex &>(side_number_table)
402 .insert(boost::shared_ptr<SideNumber>(
403 new SideNumber(ent, side_number, sense, offset)))
404 .first;
405 // std::cerr << side_number << " " << sense << " " << offset << std::endl;
406 return *miit;
407}
408std::ostream &operator<<(std::ostream &os, const RefElementFace &e) {
409 os << *e.sPtr;
410 return os;
411}
413 const boost::shared_ptr<RefEntity> &ref_ents_ptr)
414 : RefElement(ref_ents_ptr) {
415 switch (ref_ents_ptr->getEntType()) {
416 case MBEDGE:
417 break;
418 default:
419 THROW_MESSAGE("this work only for TRIs");
420 }
421}
422const boost::shared_ptr<SideNumber> &
424 moab::Interface &moab = getRefEntityPtr()->getBasicDataPtr()->moab;
425 SideNumber_multiIndex::iterator miit = side_number_table.find(ent);
426 if (miit != side_number_table.end())
427 return *miit;
428 if (sPtr->ent == ent) {
429 miit =
431 .insert(boost::shared_ptr<SideNumber>(new SideNumber(ent, 0, 0, 0)))
432 .first;
433 return *miit;
434 }
435 if (type_from_handle(ent) == MBENTITYSET) {
436 miit =
438 .insert(boost::shared_ptr<SideNumber>(new SideNumber(ent, 0, 0, 0)))
439 .first;
440 return *miit;
441 }
442 int side_number, sense, offset;
443 rval = moab.side_number(sPtr->ent, ent, side_number, sense, offset);
445 miit = const_cast<SideNumber_multiIndex &>(side_number_table)
446 .insert(boost::shared_ptr<SideNumber>(
447 new SideNumber(ent, side_number, sense, offset)))
448 .first;
449 // std::cerr << side_number << " " << sense << " " << offset << std::endl;
450 return *miit;
451}
452std::ostream &operator<<(std::ostream &os, const RefElement_EDGE &e) {
453 os << *e.sPtr;
454 return os;
455}
457 const boost::shared_ptr<RefEntity> &ref_ents_ptr)
458 : RefElement(ref_ents_ptr) {
459 switch (ref_ents_ptr->getEntType()) {
460 case MBVERTEX:
461 break;
462 default:
463 THROW_MESSAGE("this works only for TRIs");
464 }
465}
466const boost::shared_ptr<SideNumber> &
468 SideNumber_multiIndex::iterator miit = side_number_table.find(ent);
469 if (miit != side_number_table.end())
470 return *miit;
471 if (sPtr->ent == ent) {
472 miit =
474 .insert(boost::shared_ptr<SideNumber>(new SideNumber(ent, 0, 0, 0)))
475 .first;
476 return *miit;
477 }
478 if (type_from_handle(ent) == MBENTITYSET) {
479 miit =
481 .insert(boost::shared_ptr<SideNumber>(new SideNumber(ent, 0, 0, 0)))
482 .first;
483 return *miit;
484 }
485 THROW_MESSAGE("no side entity for vertex if its is not an vertex itself");
486 return nullSideNumber;
487}
488std::ostream &operator<<(std::ostream &os, const RefElement_VERTEX &e) {
489 os << *e.sPtr;
490 return os;
491}
492
493} // namespace MoFEM
#define NOT_USED(x)
Definition: definitions.h:242
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:541
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:561
multi_index_container< boost::shared_ptr< SideNumber >, indexed_by< ordered_unique< member< SideNumber, EntityHandle, &SideNumber::ent > >, ordered_non_unique< composite_key< SideNumber, const_mem_fun< SideNumber, EntityType, &SideNumber::getEntType >, member< SideNumber, signed char, &SideNumber::side_number > > > > > SideNumber_multiIndex
SideNumber_multiIndex for SideNumber.
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
std::bitset< BITREFEDGES_SIZE > BitRefEdges
Definition: Types.hpp:34
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition: Templates.hpp:1869
std::ostream & operator<<(std::ostream &os, const EntitiesFieldData::EntData &e)
keeps data about abstract EDGE finite element
RefElement_EDGE(const boost::shared_ptr< RefEntity > &ref_ents_ptr)
RefElement_MESHSET(const boost::shared_ptr< RefEntity > &ref_ents_ptr)
RefElement_PRISM(const boost::shared_ptr< RefEntity > &ref_ents_ptr)
keeps data about abstract VERTEX finite element
RefElement_VERTEX(const boost::shared_ptr< RefEntity > &ref_ents_ptr)
keeps data about abstract TRI finite element
RefElementFace(const boost::shared_ptr< RefEntity > &ref_ents_ptr)
keeps data about abstract refined finite element
virtual const BitRefEdges & getBitRefEdges() const
RefElement(const boost::shared_ptr< RefEntity > &ref_ents_ptr)
static BitRefEdges DummyBitRefEdges
boost::shared_ptr< RefEntity > & getRefEntityPtr() const
Get pointer to RefEntity.
SideNumber_multiIndex side_number_table
static const boost::shared_ptr< SideNumber > nullSideNumber
keeps data about abstract TET finite element
RefElementVolume(const boost::shared_ptr< RefEntity > &ref_ents_ptr)
const BitRefEdges & getBitRefEdges() const
Struct keeps handle to refined handle.
keeps information about side number for the finite element
boost::shared_ptr< SideNumber > getSideNumberPtr() const