v0.13.1
FieldMultiIndices.cpp
Go to the documentation of this file.
1/** \file FieldMultiIndices.cpp
2 * \brief Multi-index containers for fields
3 */
4
5
6namespace MoFEM {
7
8// Not partitioned
9const bool Idx_mi_tag::IamNotPartitioned = true;
10
11// This tag is used for partitioned problems
14
15// fields
16Field::Field(moab::Interface &moab, const EntityHandle meshset)
17 : moab(moab), meshSet(meshset), tagId(NULL), tagSpaceData(NULL),
18 tagNbCoeffData(NULL), tagName(NULL), tagNameSize(0) {
19
20 auto get_tag_data_ptr = [&](const auto name, auto &tag_data) {
22 Tag th;
23 CHKERR moab.tag_get_handle(name, th);
24 CHKERR moab.tag_get_by_ptr(th, &meshset, 1, (const void **)&tag_data);
26 };
27
28 // id
29 ierr = get_tag_data_ptr("_FieldId", tagId);
30 CHKERRABORT(PETSC_COMM_SELF, ierr);
31 // space
32 ierr = get_tag_data_ptr("_FieldSpace", tagSpaceData);
33 CHKERRABORT(PETSC_COMM_SELF, ierr);
34
35 // approx. base
36 ierr = get_tag_data_ptr("_FieldBase", tagBaseData);
37 CHKERRABORT(PETSC_COMM_SELF, ierr);
38
39 // name
40 Tag th_field_name;
41 CHKERR moab.tag_get_handle("_FieldName", th_field_name);
42 CHKERR moab.tag_get_by_ptr(th_field_name, &meshSet, 1,
43 (const void **)&tagName, &tagNameSize);
44 // name prefix
45 Tag th_field_name_data_name_prefix;
46 CHKERR moab.tag_get_handle("_FieldName_DataNamePrefix",
47 th_field_name_data_name_prefix);
48 CHKERR moab.tag_get_by_ptr(th_field_name_data_name_prefix, &meshSet, 1,
49 (const void **)&tagNamePrefixData,
51 std::string name_data_prefix((char *)tagNamePrefixData, tagNamePrefixSize);
52
53 // rank
54 std::string Tag_rank_name = "_Field_Rank_" + getName();
55 CHKERR moab.tag_get_handle(Tag_rank_name.c_str(), th_FieldRank);
56 CHKERR moab.tag_get_by_ptr(th_FieldRank, &meshSet, 1,
57 (const void **)&tagNbCoeffData);
58
59 auto get_all_tags = [&]() {
61 // order
62 ApproximationOrder def_approx_order = -1;
63 std::string tag_approximation_order_name = "_App_Order_" + getName();
64 rval = moab.tag_get_handle(tag_approximation_order_name.c_str(), 1,
65 MB_TYPE_INTEGER, th_AppOrder,
66 MB_TAG_CREAT | MB_TAG_SPARSE, &def_approx_order);
67 if (rval == MB_ALREADY_ALLOCATED)
68 rval = MB_SUCCESS;
70
71 // data
72 std::string tag_data_name = name_data_prefix + getName();
73 const int def_len = 0;
74 rval = moab.tag_get_handle(
75 tag_data_name.c_str(), def_len, MB_TYPE_DOUBLE, th_FieldData,
76 MB_TAG_CREAT | MB_TAG_VARLEN | MB_TAG_SPARSE, NULL);
77 if (rval == MB_ALREADY_ALLOCATED)
78 rval = MB_SUCCESS;
80
81 std::string tag_data_name_verts = name_data_prefix + getName() + "_V";
82 rval = moab.tag_get_handle(tag_data_name_verts.c_str(), th_FieldDataVerts);
83 if (rval == MB_SUCCESS)
85 else {
86 // Since vertex tag is not it mesh that tag is not dense, it is sparse,
87 // sinc it is set to all vertices on the mesh. Is unlikely that mesh has
88 // no vertices, then above assumption does not hold.
89 tagFieldDataVertsType = MB_TAG_SPARSE;
90 VectorDouble def_vert_data(*tagNbCoeffData);
91 def_vert_data.clear();
92 rval = moab.tag_get_handle(tag_data_name_verts.c_str(), *tagNbCoeffData,
93 MB_TYPE_DOUBLE, th_FieldDataVerts,
94 MB_TAG_CREAT | tagFieldDataVertsType,
95 &*def_vert_data.begin());
96 if (rval == MB_ALREADY_ALLOCATED)
97 rval = MB_SUCCESS;
99 }
100
102 };
103
104 auto get_all_tags_deprecated = [&]() {
106 // order
107 ApproximationOrder def_approx_order = -1;
108 std::string tag_approximation_order_name = "_App_Order_" + getName();
109 rval = moab.tag_get_handle(tag_approximation_order_name.c_str(), 1,
110 MB_TYPE_INTEGER, th_AppOrder,
111 MB_TAG_CREAT | MB_TAG_SPARSE, &def_approx_order);
112 if (rval == MB_ALREADY_ALLOCATED)
113 rval = MB_SUCCESS;
115
116 // data
117 std::string tag_data_name = name_data_prefix + getName();
118 const int def_len = 0;
119 rval = moab.tag_get_handle(
120 tag_data_name.c_str(), def_len, MB_TYPE_DOUBLE, th_FieldData,
121 MB_TAG_CREAT | MB_TAG_VARLEN | MB_TAG_SPARSE, NULL);
122 if (rval == MB_ALREADY_ALLOCATED)
123 rval = MB_SUCCESS;
125
126 std::string tag_data_name_verts = name_data_prefix + getName() + "V";
127 rval = moab.tag_get_handle(tag_data_name_verts.c_str(), th_FieldDataVerts);
128 if (rval == MB_SUCCESS)
130 else {
131 // Since vertex tag is not it mesh that tag is not dense, it is sparse,
132 // sinc it is set to all vertices on the mesh. Is unlikely that mesh has
133 // no vertices, then above assumption does not hold.
134 tagFieldDataVertsType = MB_TAG_SPARSE;
135 VectorDouble def_vert_data(*tagNbCoeffData);
136 def_vert_data.clear();
137 rval = moab.tag_get_handle(tag_data_name_verts.c_str(), *tagNbCoeffData,
138 MB_TYPE_DOUBLE, th_FieldDataVerts,
139 MB_TAG_CREAT | tagFieldDataVertsType,
140 &*def_vert_data.begin());
141 if (rval == MB_ALREADY_ALLOCATED)
142 rval = MB_SUCCESS;
144 }
145
147 };
148
149 Version file_ver;
151 CHK_THROW_MESSAGE(ierr, "Not known file version");
152 if (file_ver.majorVersion >= 0 && file_ver.minorVersion >= 12 &&
153 file_ver.buildVersion >= 1) {
154 ierr = get_all_tags();
155 CHKERRABORT(PETSC_COMM_SELF, ierr);
156 } else {
157 ierr = get_all_tags_deprecated();
158 CHKERRABORT(PETSC_COMM_SELF, ierr);
159 }
160
162
163 auto reset_entity_order_table = [&]() {
164 for (int tt = 0; tt != MBMAXTYPE; ++tt)
165 forderTable[tt] = NULL;
166 };
167
168 auto set_entity_order_table = [&]() {
169 switch (*tagBaseData) {
172 switch (*tagSpaceData) {
173 case H1:
174 forderTable[MBVERTEX] = [](int P) -> int { return (P > 0) ? 1 : 0; };
175 forderTable[MBEDGE] = [](int P) -> int { return NBEDGE_H1(P); };
176 forderTable[MBTRI] = [](int P) -> int { return NBFACETRI_H1(P); };
177 forderTable[MBQUAD] = [](int P) -> int { return NBFACEQUAD_H1(P); };
178 forderTable[MBTET] = [](int P) -> int { return NBVOLUMETET_H1(P); };
179 forderTable[MBHEX] = [](int P) -> int { return NBVOLUMEHEX_H1(P); };
180 forderTable[MBPRISM] = [](int P) -> int { return NBVOLUMEPRISM_H1(P); };
181 break;
182 case HCURL:
183 forderTable[MBVERTEX] = [](int P) -> int {
184 (void)P;
185 return 0;
186 };
187 forderTable[MBEDGE] = [](int P) -> int {
188 return NBEDGE_AINSWORTH_HCURL(P);
189 };
190 forderTable[MBTRI] = [](int P) -> int {
192 };
193 forderTable[MBTET] = [](int P) -> int {
195 };
196 break;
197 case HDIV:
198 forderTable[MBVERTEX] = [](int P) -> int {
199 (void)P;
200 return 0;
201 };
202 forderTable[MBEDGE] = [](int P) -> int {
203 (void)P;
204 return NBEDGE_HDIV(P);
205 };
206 forderTable[MBTRI] = [](int P) -> int {
207 return NBFACETRI_AINSWORTH_HDIV(P);
208 };
209 forderTable[MBTET] = [](int P) -> int {
211 };
212 break;
213 case L2:
214 forderTable[MBVERTEX] = [](int P) -> int {
215 (void)P;
216 return 1;
217 };
218 forderTable[MBEDGE] = [](int P) -> int { return NBEDGE_L2(P); };
219 forderTable[MBTRI] = [](int P) -> int { return NBFACETRI_L2(P); };
220 forderTable[MBQUAD] = [](int P) -> int { return NBFACEQUAD_L2(P); };
221 forderTable[MBTET] = [](int P) -> int { return NBVOLUMETET_L2(P); };
222 forderTable[MBHEX] = [](int P) -> int { return NBVOLUMEHEX_L2(P); };
223 break;
224 default:
225 THROW_MESSAGE("unknown approximation space");
226 }
227 break;
229 switch (*tagSpaceData) {
230 case H1:
231 forderTable[MBVERTEX] = [](int P) -> int { return (P > 0) ? 1 : 0; };
232 forderTable[MBEDGE] = [](int P) -> int { return NBEDGE_H1(P); };
233 forderTable[MBTRI] = [](int P) -> int { return NBFACETRI_H1(P); };
234 forderTable[MBQUAD] = [](int P) -> int { return NBFACEQUAD_H1(P); };
235 forderTable[MBTET] = [](int P) -> int { return NBVOLUMETET_H1(P); };
236 forderTable[MBPRISM] = [](int P) -> int { return NBVOLUMEPRISM_H1(P); };
237 break;
238 case L2:
239 forderTable[MBVERTEX] = [](int P) -> int {
240 (void)P;
241 return 1;
242 };
243 forderTable[MBEDGE] = [](int P) -> int { return NBEDGE_L2(P); };
244 forderTable[MBTRI] = [](int P) -> int { return NBFACETRI_L2(P); };
245 forderTable[MBQUAD] = [](int P) -> int { return NBFACEQUAD_L2(P); };
246 forderTable[MBTET] = [](int P) -> int { return NBVOLUMETET_L2(P); };
247 forderTable[MBHEX] = [](int P) -> int { return NBVOLUMEHEX_L2(P); };
248 break;
249 default:
250 THROW_MESSAGE("unknown approximation space or not yet implemented");
251 }
252 break;
254 switch (*tagSpaceData) {
255 case H1:
256 forderTable[MBVERTEX] = [](int P) -> int { return (P > 0) ? 1 : 0; };
257 forderTable[MBEDGE] = [](int P) -> int { return NBEDGE_H1(P); };
258 forderTable[MBTRI] = [](int P) -> int { return NBFACETRI_H1(P); };
259 forderTable[MBQUAD] = [](int P) -> int { return NBFACEQUAD_H1(P); };
260 forderTable[MBTET] = [](int P) -> int { return NBVOLUMETET_H1(P); };
261 forderTable[MBHEX] = [](int P) -> int { return NBVOLUMEHEX_H1(P); };
262 forderTable[MBPRISM] = [](int P) -> int { return NBVOLUMEPRISM_H1(P); };
263 break;
264 case HCURL:
265 forderTable[MBVERTEX] = [](int P) -> int {
266 (void)P;
267 return 0;
268 };
269 forderTable[MBEDGE] = [](int P) -> int {
270 return NBEDGE_DEMKOWICZ_HCURL(P);
271 };
272 forderTable[MBTRI] = [](int P) -> int {
274 };
275 forderTable[MBQUAD] = [](int P) -> int {
277 };
278 forderTable[MBTET] = [](int P) -> int {
280 };
281 forderTable[MBHEX] = [](int P) -> int {
283 };
284 break;
285 case HDIV:
286 forderTable[MBVERTEX] = [](int P) -> int {
287 (void)P;
288 return 0;
289 };
290 forderTable[MBEDGE] = [](int P) -> int {
291 (void)P;
292 return 0;
293 };
294 forderTable[MBTRI] = [](int P) -> int {
295 return NBFACETRI_DEMKOWICZ_HDIV(P);
296 };
297 forderTable[MBQUAD] = [](int P) -> int {
299 };
300 forderTable[MBTET] = [](int P) -> int {
302 };
303 forderTable[MBHEX] = [](int P) -> int {
305 };
306 break;
307 case L2:
308 forderTable[MBVERTEX] = [](int P) -> int {
309 (void)P;
310 return 1;
311 };
312 forderTable[MBEDGE] = [](int P) -> int { return NBEDGE_L2(P); };
313 forderTable[MBTRI] = [](int P) -> int { return NBFACETRI_L2(P); };
314 forderTable[MBQUAD] = [](int P) -> int { return NBFACEQUAD_L2(P); };
315 forderTable[MBTET] = [](int P) -> int { return NBVOLUMETET_L2(P); };
316 forderTable[MBHEX] = [](int P) -> int { return NBVOLUMEHEX_L2(P); };
317 break;
318 default:
319 THROW_MESSAGE("unknown approximation space or not yet implemented");
320 }
321 break;
322 case USER_BASE:
323 for (int ee = 0; ee < MBMAXTYPE; ee++) {
324 forderTable[ee] = [](int P) -> int {
325 (void)P;
326 return 0;
327 };
328 }
329 break;
330 default:
331 if (*tagSpaceData != NOFIELD) {
332 THROW_MESSAGE("unknown approximation base");
333 } else {
334 for (EntityType t = MBVERTEX; t < MBMAXTYPE; t++)
335 forderTable[t] = [](int P) -> int {
336 (void)P;
337 return 1;
338 };
339 }
340 }
341 };
342
343 reset_entity_order_table();
344 set_entity_order_table();
346 CHKERRABORT(PETSC_COMM_SELF, ierr);
347};
348
351
352 for (auto t = MBVERTEX; t != MBMAXTYPE; ++t) {
353
354 int DD = 0;
355 int nb_last_order_dofs = 0;
356 const int rank = (*tagNbCoeffData);
357 if (forderTable[t]) {
358
359 for (int oo = 0; oo < MAX_DOFS_ON_ENTITY; ++oo) {
360
361 const int nb_order_dofs = forderTable[t](oo);
362 const int diff_oo = nb_order_dofs - nb_last_order_dofs;
363 if (diff_oo >= 0) {
364
365 if ((DD + rank * diff_oo) < MAX_DOFS_ON_ENTITY)
366 for (int dd = 0; dd < diff_oo; ++dd)
367 for (int rr = 0; rr != rank; ++rr, ++DD)
368 dofOrderMap[t][DD] = oo;
369 else
370 break;
371
372 nb_last_order_dofs = nb_order_dofs;
373
374 } else {
375 break;
376 }
377 }
378 }
379
380 std::fill(&dofOrderMap[t][DD], dofOrderMap[t].end(), -1);
381 }
382
384}
385
386std::ostream &operator<<(std::ostream &os, const Field &e) {
387 os << e.getNameRef() << " field_id " << e.getId().to_ulong() << " space "
388 << FieldSpaceNames[e.getSpace()] << " approximation base "
389 << ApproximationBaseNames[e.getApproxBase()] << " nb coefficients "
390 << e.getNbOfCoeffs() << " meshset " << e.meshSet;
391 return os;
392}
393
394// FieldEntityEntFiniteElementAdjacencyMap
397 const boost::shared_ptr<FieldEntity> &ent_field_ptr,
398 const boost::shared_ptr<EntFiniteElement> &ent_fe_ptr)
399 : byWhat(0), entFieldPtr(ent_field_ptr), entFePtr(ent_fe_ptr) {}
400
401std::ostream &operator<<(std::ostream &os,
403 os << "byWhat " << std::bitset<3>(e.byWhat) << " " << *e.entFieldPtr
404 << std::endl
405 << *e.entFePtr->getFiniteElementPtr();
406 return os;
407}
408
409} // namespace MoFEM
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:236
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:62
@ USER_BASE
user implemented approximation base
Definition: definitions.h:68
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:541
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:595
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ NOFIELD
scalar or vector of scalars describe (no true field)
Definition: definitions.h:84
@ H1
continuous field
Definition: definitions.h:85
@ HCURL
field with continuous tangents
Definition: definitions.h:86
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
static const char *const FieldSpaceNames[]
Definition: definitions.h:92
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
static const char *const ApproximationBaseNames[]
Definition: definitions.h:72
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:561
#define NBFACEQUAD_DEMKOWICZ_HDIV(P)
#define NBEDGE_DEMKOWICZ_HCURL(P)
#define NBVOLUMETET_H1(P)
Number of base functions on tetrahedron for H1 space.
#define NBVOLUMETET_AINSWORTH_HDIV(P)
#define NBVOLUMEHEX_DEMKOWICZ_HCURL(P)
#define NBFACEQUAD_H1(P)
Number of base functions on quad for H1 space.
#define NBVOLUMETET_AINSWORTH_HCURL(P)
#define NBFACETRI_AINSWORTH_HCURL(P)
#define NBFACEQUAD_L2(P)
Number of base functions on quad for L2 space.
#define NBVOLUMEHEX_H1(P)
Number of base functions on hex for H1 space.
#define NBFACETRI_L2(P)
Number of base functions on triangle for L2 space.
#define NBVOLUMETET_DEMKOWICZ_HDIV(P)
#define NBEDGE_HDIV(P)
#define NBVOLUMEPRISM_H1(P)
Number of base functions on prism for H1 space.
#define NBVOLUMEHEX_DEMKOWICZ_HDIV(P)
#define NBFACEQUAD_DEMKOWICZ_HCURL(P)
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.
#define NBFACETRI_DEMKOWICZ_HDIV(P)
#define NBEDGE_L2(P)
Number of base functions on edge fro L2 space.
#define NBVOLUMEHEX_L2(P)
Number of base functions on hexahedron for L2 space.
#define NBVOLUMETET_DEMKOWICZ_HCURL(P)
#define NBFACETRI_DEMKOWICZ_HCURL(P)
#define NBEDGE_AINSWORTH_HCURL(P)
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.
#define NBFACETRI_AINSWORTH_HDIV(P)
#define NBVOLUMETET_L2(P)
Number of base functions on tetrahedron for L2 space.
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
int ApproximationOrder
Approximation on the entity.
Definition: Types.hpp:26
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
std::ostream & operator<<(std::ostream &os, const EntitiesFieldData::EntData &e)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
constexpr double t
plate stiffness
Definition: plate.cpp:59
FieldEntityEntFiniteElementAdjacencyMap of mofem finite element and entities.
const boost::shared_ptr< FieldEntity > entFieldPtr
field entity
FieldEntityEntFiniteElementAdjacencyMap(const boost::shared_ptr< FieldEntity > &ent_field_ptr, const boost::shared_ptr< EntFiniteElement > &ent_fe_ptr)
const boost::shared_ptr< EntFiniteElement > entFePtr
finite element entity
Provide data structure for (tensor) field approximation.
FieldCoefficientsNumber getNbOfCoeffs() const
Get number of field coefficients.
const void * tagNamePrefixData
tag keeps name prefix of the field
FieldOrderTable forderTable
nb. DOFs table for entities
const void * tagName
tag keeps name of the field
int tagNameSize
number of bits necessary to keep field name
std::string getName() const
Get field name.
DofsOrderMap dofOrderMap
FieldApproximationBase getApproxBase() const
Get approximation base.
unsigned int bitNumber
EntityHandle meshSet
keeps entities for this meshset
FieldSpace * tagSpaceData
tag keeps field space
FieldBitNumber getBitNumberCalculate() const
Calculate number of set bit in Field ID. Each field has uid, get getBitNumber get number of bit set f...
Tag th_FieldData
Tag storing field values on entity in the field.
MoFEMErrorCode rebuildDofsOrderMap()
FieldCoefficientsNumber * tagNbCoeffData
BitFieldId * tagId
Tag field rank.
FieldSpace getSpace() const
Get field approximation space.
Field(moab::Interface &moab, const EntityHandle meshset)
constructor for moab field
Tag th_FieldDataVerts
Tag storing field values on vertices in the field.
FieldApproximationBase * tagBaseData
tag keeps field base
TagType tagFieldDataVertsType
Tag th_AppOrder
Tag storing approximation order on entity.
moab::Interface & moab
boost::string_ref getNameRef() const
Get string reference to field name.
const BitFieldId & getId() const
Get unique field id.
static const bool IamNotPartitioned
static const bool IamNotPartitioned
static const bool IamNotPartitioned
static MoFEMErrorCode getFileVersion(moab::Interface &moab, Version &version)
Get database major version.