88 {
89 int count_meshset_sides = 0;
91
92 auto check = [](
auto &
b,
auto &
m,
auto &
bit) {
93 return ((
bit & b).any()) && ((
bit &
m) ==
bit);
94 };
95
96 auto set_child_data_entity = [](auto &parent_data, auto &child_data) {
98 child_data.getEntDataBitRefLevel() = parent_data.getEntDataBitRefLevel();
99 child_data.sPace = parent_data.getSpace();
100 child_data.bAse = parent_data.getBase();
101 child_data.sEnse = parent_data.getSense();
102 child_data.oRder = parent_data.getOrder();
103 child_data.iNdices.swap(parent_data.getIndices());
104 child_data.localIndices.swap(parent_data.getLocalIndices());
105 child_data.dOfs.swap(parent_data.getFieldDofs());
106 child_data.fieldData.swap(parent_data.getFieldData());
107 child_data.fieldEntities.swap(parent_data.getFieldEntities());
109 };
110
111 auto set_child_data_vertex = [&](auto &parent_data, auto &child_data,
112 int node, int num_nodes) {
114
115 child_data.getEntDataBitRefLevel().resize(1, false);
116 child_data.getEntDataBitRefLevel()[0] =
117 parent_data.getEntDataBitRefLevel()[node];
118 child_data.sPace = parent_data.getSpace();
119 child_data.bAse = parent_data.getBase();
120 child_data.sEnse = parent_data.getSense();
121 child_data.oRder = parent_data.getOrder();
122
123 if (parent_data.dOfs.size()) {
124
125 auto &field_entities = parent_data.getFieldEntities();
126 auto &vertex_entity = field_entities[node];
127
128 const auto nb_coeffs = vertex_entity->getNbOfCoeffs();
129
130 child_data.fieldEntities.resize(1, false);
131 child_data.fieldEntities[0] = vertex_entity;
132#ifndef NDEBUG
133 if (parent_data.dOfs.size() != nb_coeffs * num_nodes)
135 "Inconsistent number of DOFs and vertices %d != %d",
136 parent_data.dOfs.size(), nb_coeffs * num_nodes);
137#endif
138
139
140
141 if (parent_data.iNdices.size()) {
142
143#ifndef NDEBUG
144 if (parent_data.dOfs.size() != parent_data.iNdices.size())
146 "Inconsistent size of DOFs %d != %d",
147 parent_data.dOfs.size(), parent_data.iNdices.size());
148#endif
149
150 child_data.iNdices.resize(nb_coeffs, false);
151 child_data.localIndices.resize(nb_coeffs, false);
152 child_data.dOfs.resize(nb_coeffs, false);
153 child_data.fieldData.resize(nb_coeffs, false);
154
155 int DD = 0;
156 for (
auto dd = node * nb_coeffs;
dd != (node + 1) * nb_coeffs;
158 child_data.iNdices[DD] = parent_data.iNdices[
dd];
159 child_data.localIndices[DD] = parent_data.localIndices[
dd];
160 child_data.dOfs[DD] = parent_data.dOfs[
dd];
161 child_data.fieldData[DD] = parent_data.fieldData[
dd];
162 }
163 } else {
164
165 child_data.iNdices.clear();
166 child_data.localIndices.clear();
167 child_data.dOfs.resize(nb_coeffs, false);
168 child_data.fieldData.resize(nb_coeffs, false);
169 int DD = 0;
170 for (
auto dd = node * nb_coeffs;
dd != (node + 1) * nb_coeffs;
172 child_data.dOfs[DD] = parent_data.dOfs[
dd];
173 child_data.fieldData[DD] = parent_data.fieldData[
dd];
174 }
175
176 }
177 }
178
180 };
181
182
183
184
185
186
187
188 auto set_child_base_entity = [](auto &parent_data, auto &child_data) {
190 child_data.resetFieldDependentData();
191 child_data.getEntDataBitRefLevel() = parent_data.getEntDataBitRefLevel();
192 child_data.sPace = parent_data.getSpace();
193
194#ifndef NDEBUG
197#endif
198
201 for (auto derivative = 0; derivative != BaseDerivatives::LastDerivative;
202 ++derivative) {
203 auto parent_base =
206 if (parent_base) {
207 auto &child_base = child_data.getNSharedPtr(
210 if (!child_base)
211 child_base = boost::make_shared<MatrixDouble>();
212 child_base->swap(*parent_base);
213 }
214 }
215 }
216
218 };
219
220 auto set_child_base_vertex = [](auto &parent_data, auto &child_data, int node,
221 int num_nodes) {
223 child_data.resetFieldDependentData();
224 child_data.getEntDataBitRefLevel().resize(1, false);
225 child_data.getEntDataBitRefLevel()[0] =
226 parent_data.getEntDataBitRefLevel()[node];
227 child_data.sPace = parent_data.getSpace();
228
229#ifndef NDEBUG
232#endif
233
236 for (auto derivative = 0; derivative != BaseDerivatives::LastDerivative;
237 ++derivative) {
238 auto parent_base =
241 if (parent_base) {
242 auto &child_base = child_data.getNSharedPtr(
245
246 const auto num_bases_per_node = parent_base->size2() / num_nodes;
247 if (parent_base->size2() % num_nodes) {
249 "Inconsistent nb of base functions and nodes mod(%d, %d)",
250 parent_base->size2(), num_nodes);
251 }
252
253 if (!child_base)
254 child_base = boost::make_shared<MatrixDouble>(parent_base->size1(),
255 num_bases_per_node);
256 else
257 child_base->resize(parent_base->size1(), num_bases_per_node, false);
258
259 for (auto gg = 0; gg != parent_base->size1(); ++gg) {
260 int DD = 0;
261 for (auto dd = node * num_bases_per_node;
262 dd != (node + 1) * num_bases_per_node; ++
dd, ++DD) {
263 (*child_base)(gg, DD) = (*parent_base)(gg,
dd);
264 }
265 }
266 }
267 }
268 }
269
271 };
272
273
274
275
276
277 auto switch_off_dofs_children = [&](auto &parent_ent_data, auto &child_data) {
279 for (
auto i : parent_ent_data.getIndices()) {
281 for (auto &child_ent_data : child_data) {
282 auto it = std::find(child_ent_data.getIndices().begin(),
283 child_ent_data.getIndices().end(),
i);
284 if (it != child_ent_data.getIndices().end()) {
285 const auto dof_idx =
286 std::distance(child_ent_data.getIndices().begin(), it);
287 auto &bit_dof =
288 child_ent_data.getFieldDofs()[dof_idx]->getBitRefLevel();
290 child_ent_data.getIndices()[dof_idx] = -1;
291 child_ent_data.getLocalIndices()[dof_idx] = -1;
292 child_ent_data.getFieldData()[dof_idx] = 0;
293 }
294 }
295 }
296 }
297 }
299 };
300
301
302
303 auto do_work_parent_hook = [&](
DataOperator *op_ptr,
int side,
305 EntitiesFieldData::EntData &data) {
307
309 auto &field_entities = data.getFieldEntities();
310
311#ifndef NDEBUG
314 << "Add entity data to meshset "
315 << "side/type: " << side << "/" << CN::EntityTypeName(type)
319 }
321 for (auto field_ent : field_entities) {
323 << "Parent field entity bit: " << field_ent->getBitRefLevel() << " "
324 << *field_ent;
325 }
326 }
327#endif
328
329 const auto nb_ents = field_entities.size();
330 const auto nb_nodes = data.getEntDataBitRefLevel().size();
331
332
333 if (type == MBVERTEX &&
334 nb_nodes != up_op_ptr->getNumberOfNodesOnElement()) {
335 SETERRQ2(
337 "Inconsistency between bit ref levels and number of nodes %d != %d",
338 nb_nodes, up_op_ptr->getNumberOfNodesOnElement());
339 }
340
341
342 auto get_child_meshset =
344 auto &data_on_meshset = entities_field_data.dataOnEntities[MBENTITYSET];
345 if (data_on_meshset.size() < count_meshset_sides) {
347
348 data_on_meshset.transfer(data_on_meshset.end(),
350 } else {
351
352 entities_field_data.dataOnEntities[MBENTITYSET].resize(
353 count_meshset_sides);
354 }
355 }
356 return entities_field_data
357 .dataOnEntities[MBENTITYSET][count_meshset_sides - 1];
358 };
359
361 for (auto node = 0; node != nb_nodes; ++node) {
364 ++count_meshset_sides;
365 auto &child_data_meshset = get_child_meshset(count_meshset_sides);
366 if (type == MBVERTEX) {
367 const auto num_nodes = up_op_ptr->getNumberOfNodesOnElement();
368 CHKERR set_child_base_vertex(data, child_data_meshset, node,
369 num_nodes);
370 } else {
371 CHKERR set_child_base_entity(data, child_data_meshset);
372 }
373 }
374 }
375
376 } else {
377
378 {
379 boost::container::static_vector<EntityType, 128> ents_type;
380 ents_type.reserve(field_entities.size());
381 for (auto &field_entity : field_entities)
382 ents_type.push_back(field_entity->getEntType());
383 std::sort(ents_type.begin(), ents_type.end());
384 auto end = std::unique(ents_type.begin(), ents_type.end());
385 for (auto it_t = ents_type.begin(); it_t != end; ++it_t)
386 CHKERR switch_off_dofs_children(
387 data, entities_field_data.dataOnEntities[*it_t]);
388 if (type == MBENTITYSET)
389 CHKERR switch_off_dofs_children(
390 data, entities_field_data.dataOnEntities[type]);
391 }
392
393 if (type == MBVERTEX &&
394 nb_ents != up_op_ptr->getNumberOfNodesOnElement()) {
395 SETERRQ2(
397 "Number of nodes is expected to much number of entities: %d != %d",
398 nb_ents, up_op_ptr->getNumberOfNodesOnElement());
399 } else if (type != MBVERTEX && nb_ents > 1) {
401 "Code can handle only one entity");
402 }
403
404 for (auto node = 0; node != nb_nodes; ++node) {
405 auto &bit_ent = data.getEntDataBitRefLevel()[node];
407 ++count_meshset_sides;
408 auto &child_data_meshset = get_child_meshset(count_meshset_sides);
409 if (type == MBVERTEX) {
410 const auto num_nodes = up_op_ptr->getNumberOfNodesOnElement();
411 CHKERR set_child_data_vertex(data, child_data_meshset, node,
412 num_nodes);
413 } else {
414 CHKERR set_child_data_entity(data, child_data_meshset);
415 }
416 }
417 }
418 }
419
421 };
422
423
426
429 }
430
431 auto loop_parent_fe = [&]() {
433
434
435 parentElePtr->getOpPtrVector().back().doWorkRhsHook = do_work_parent_hook;
438
439#ifndef NDEBUG
442 if (
getGaussPts().size1() != parent_gauss_pts.size1()) {
444 << "Calling element: "
445 << boost::typeindex::type_id_runtime(*parentElePtr).pretty_name();
447 MOFEM_LOG(
"SELF", Sev::error) << parent_gauss_pts;
449 "Wrong number of weights");
450 }
451 if (
getGaussPts().size2() != parent_gauss_pts.size2()) {
453 << "Calling element: "
454 << boost::typeindex::type_id_runtime(*parentElePtr).pretty_name();
456 "Wrong number of integration points");
457 }
458 }
459#endif
460
462 };
463
464
466 }
467
470 << "Number of meshset entities "
471 << entities_field_data.dataOnEntities[MBENTITYSET].size();
472 }
473
474 auto &data_on_meshset = entities_field_data.dataOnEntities[MBENTITYSET];
475 auto it = data_on_meshset.begin();
476
477
478 if (count_meshset_sides < data_on_meshset.size()) {
479 for (auto s = 0; s != count_meshset_sides; ++s)
480 ++it;
482 data_on_meshset);
483 }
484
485 auto set_up_derivative_ent_data = [&](auto &entities_field_data,
486 auto &derivative_entities_field_data) {
488
489 using EntData = EntitiesFieldData::EntData;
490 using DerivedEntData = DerivedEntitiesFieldData::DerivedEntData;
491
493 auto &ents_data = entities_field_data.dataOnEntities[MBENTITYSET];
494 auto &dev_ents_data =
495 derivative_entities_field_data.dataOnEntities[MBENTITYSET];
496 dev_ents_data.clear();
497 for (
auto c = 0;
c < ents_data.size(); ++
c) {
498 boost::shared_ptr<EntData> ent_data_ptr(data_ptr, &ents_data[
c]);
499 dev_ents_data.push_back(new DerivedEntData(ent_data_ptr));
500 }
502 };
503
505 CHKERR set_up_derivative_ent_data(
506 entities_field_data,
508 }
509
510#ifndef NDEBUG
511 auto &side_table =
513 for (auto &data : entities_field_data.dataOnEntities[MBENTITYSET]) {
514 for (auto dof : data.getFieldDofs()) {
515 auto &bit_dof = dof->getBitRefLevel();
517 auto ent = dof->getEnt();
518 if (side_table.find(ent) == side_table.end()) {
520 "Adjacency not found");
521 }
522 }
523 }
524 }
525#endif
526
528}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
static const char *const FieldSpaceNames[]
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
static const char *const ApproximationBaseNames[]
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
FTensor::Index< 'm', SPACE_DIM > m
#define MOFEM_LOG(channel, severity)
Log.
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
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)
EntitiesFieldData::EntData::BaseDerivatives BaseDerivatives
DataOperator(const bool symm=true)
Data on single entity (This is passed as argument to DataOperator::doWork)
virtual std::vector< BitRefLevel > & getEntDataBitRefLevel()
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
MoFEMErrorCode loopParent(const string &fe_name, ForcesAndSourcesCore *parent_fe, const int verb=QUIET, const LogManager::SeverityLevel sev=Sev::noisy)
User calls this function to loop over parent elements. This function calls finite element with its op...
std::string getFEName() const
Get name of the element.
ForcesAndSourcesCore * getPtrFE() const
@ OPSPACE
operator do Work is execute on space data
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
auto & getDataOnElementBySpaceArray()
Get data on entities and space.
boost::ptr_deque< EntitiesFieldData::EntData > poolEntsVector
Public Member Functions inherited from MoFEM::ForcesAndSourcesCore::UserDataOperator