v0.13.2
Loading...
Searching...
No Matches
MeshProjectionDataOperators.cpp
Go to the documentation of this file.
1/** \file MeshProjectionDataOperators.cpp
2
3
4*/
5
6
7
8namespace MoFEM {
9
10OpRunParent::OpRunParent(boost::shared_ptr<ForcesAndSourcesCore> parent_ele_ptr,
11 BitRefLevel bit_parent, BitRefLevel bit_parent_mask,
12 boost::shared_ptr<ForcesAndSourcesCore> this_ele_ptr,
13 BitRefLevel bit_this, BitRefLevel bit_this_mask,
14 int verb, Sev sev)
17 parentElePtr(parent_ele_ptr), bitParent(bit_parent),
18 bitParentMask(bit_parent_mask), thisElePtr(this_ele_ptr),
19 bitThis(bit_this), bitThisMask(bit_this_mask), verbosity(verb),
20 severityLevel(sev) {}
21
25
26 auto &bit = getFEMethod()->numeredEntFiniteElementPtr->getBitRefLevel();
27
28 auto check = [&](auto &b, auto &m) {
29 return ((bit & b).any()) && ((bit & m) == bit);
30 };
31
32#ifndef NDEBUG
33 if (verbosity > QUIET) {
34 MOFEM_LOG_CHANNEL("SELF");
35 MOFEM_TAG_AND_LOG("SELF", severityLevel, "OpRunParent")
36 << "FE bit " << bit
37 << " check parent = " << check(bitParent, bitParentMask)
38 << " check this " << check(bitThis, bitThisMask);
39 }
40#endif
41
42 if (check(bitParent, bitParentMask)) {
43 if (parentElePtr)
46
47 } else if (check(bitThis, bitThisMask)) {
48 if (thisElePtr)
50 }
51
53}
54
56 std::string field_name, OpType op_parent_type,
57 boost::shared_ptr<ForcesAndSourcesCore> parent_ele_ptr,
58 BitRefLevel bit_child, BitRefLevel bit_child_mask,
59 BitRefLevel bit_parent_ent, BitRefLevel bit_parent_ent_mask, int verb,
60 Sev sev)
62 fieldName(field_name), opParentType(op_parent_type),
63 parentElePtr(parent_ele_ptr), bitChild(bit_child),
64 bitChildMask(bit_child_mask), bitParentEnt(bit_parent_ent),
65 bitParentEntMask(bit_parent_ent_mask), verbosity(verb),
66 severityLevel(sev) {
67 // Push op to collect data
68 auto field_op =
70 parentElePtr->getOpPtrVector().push_back(field_op);
71}
72
74 FieldSpace space, OpType op_parent_type,
75 boost::shared_ptr<ForcesAndSourcesCore> parent_ele_ptr,
76 BitRefLevel bit_child, BitRefLevel bit_child_mask,
77 BitRefLevel bit_parent_ent, BitRefLevel bit_parent_ent_mask, int verb,
78 Sev sev)
79 : ForcesAndSourcesCore::UserDataOperator(space, op_parent_type),
80 fieldName(""), approxSpace(space), opParentType(op_parent_type),
81 parentElePtr(parent_ele_ptr), bitChild(bit_child),
82 bitChildMask(bit_child_mask), bitParentEnt(bit_parent_ent),
83 bitParentEntMask(bit_parent_ent_mask), verbosity(verb),
84 severityLevel(sev) {
85 // Push op to collect data
86 auto field_op =
88 parentElePtr->getOpPtrVector().push_back(field_op);
89}
90
92 const bool error_if_no_base) {
93 int count_meshset_sides = 0;
95
96 auto check = [](auto &b, auto &m, auto &bit) {
97 return ((bit & b).any()) && ((bit & m) == bit);
98 };
99
100 auto set_child_data = [](auto &parent_data, auto &child_data) {
102 child_data.getEntDataBitRefLevel() = parent_data.getEntDataBitRefLevel();
103 child_data.sPace = parent_data.getSpace();
104 child_data.bAse = parent_data.getBase();
105 child_data.sEnse = parent_data.getSense();
106 child_data.oRder = parent_data.getOrder();
107 child_data.iNdices.swap(parent_data.getIndices());
108 child_data.localIndices.swap(parent_data.getLocalIndices());
109 child_data.dOfs.swap(parent_data.getFieldDofs());
110 child_data.fieldData.swap(parent_data.getFieldData());
111 child_data.fieldEntities.swap(parent_data.getFieldEntities());
113 };
114
115 /**
116 * @brief swap child base
117 *
118 * @todo add swap for Bernstein-Bezier base
119 *
120 */
121 auto set_child_base = [](auto &parent_data, auto &child_data) {
123 child_data.resetFieldDependentData();
124 child_data.getEntDataBitRefLevel() = parent_data.getEntDataBitRefLevel();
125 child_data.sPace = parent_data.getSpace();
126 child_data.getEntDataBitRefLevel() = parent_data.getEntDataBitRefLevel();
127
128#ifndef NDEBUG
129 if (child_data.bAse == AINSWORTH_BERNSTEIN_BEZIER_BASE)
130 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Base not implemented");
131#endif
132
134 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
135 for (auto derivative = 0; derivative != BaseDerivatives::LastDerivative;
136 ++derivative) {
137 auto parent_base =
138 parent_data.getNSharedPtr(static_cast<FieldApproximationBase>(b),
139 static_cast<BaseDerivatives>(derivative));
140 if (parent_base) {
141 auto child_base = child_data.getNSharedPtr(
142 static_cast<FieldApproximationBase>(b),
143 static_cast<BaseDerivatives>(derivative));
144 if (!child_base)
145 child_base = boost::make_shared<MatrixDouble>();
146 child_base->swap(*parent_base);
147 }
148 }
149 }
150
152 };
153
154 auto switch_of_dofs_children = [](auto &parent_ent_data, auto &child_data) {
156 for (auto i : parent_ent_data.getIndices()) {
157 if (i >= 0) {
158 for (auto &child_ent_data : child_data) {
159 auto it = std::find(child_ent_data.getIndices().begin(),
160 child_ent_data.getIndices().end(), i);
161 if (it != child_ent_data.getIndices().end()) {
162 const auto dof_idx =
163 std::distance(child_ent_data.getIndices().begin(), it);
164 child_ent_data.getIndices()[dof_idx] = -1;
165 child_ent_data.getLocalIndices()[dof_idx] = -1;
166 child_ent_data.getFieldData()[dof_idx] = 0;
167 }
168 }
169 }
170 }
172 };
173
174 // that forces to run operator at last instance and collect data on
175 // entities
176 auto do_work_parent_hook = [&](DataOperator *op_ptr, int side,
177 EntityType type,
180
181 BitRefLevel &bit_ent = data.getEntDataBitRefLevel();
182
183 // note all nodes from all added
184 if (check(bitParentEnt, bitParentEntMask, bit_ent)) {
185 ++count_meshset_sides;
186
187 if (verbosity >= VERBOSE) {
188 MOFEM_LOG("SELF", severityLevel)
189 << "Add entity data to meshset "
190 << "side/type: " << side << "/" << CN::EntityTypeName(type)
191 << " op space/base: " << FieldSpaceNames[data.getSpace()] << "/"
192 << ApproximationBaseNames[data.getBase()];
193 }
194
195 auto &data_on_meshset = entities_field_data.dataOnEntities[MBENTITYSET];
196 if (data_on_meshset.size() < count_meshset_sides) {
197 if (poolEntsVector.size()) {
198 data_on_meshset.transfer(data_on_meshset.end(),
200 } else {
201 entities_field_data.dataOnEntities[MBENTITYSET].resize(
202 count_meshset_sides);
203 }
204 }
205
206 auto &child_data_meshset =
207 entities_field_data
208 .dataOnEntities[MBENTITYSET][count_meshset_sides - 1];
209
210 if (opParentType == OPSPACE) {
211 CHKERR set_child_base(data, child_data_meshset);
212 child_data_meshset.resetFieldDependentData();
213
214 } else {
215
216 auto &field_entities = data.getFieldEntities();
217
218 // note all nodes from all added
219 if (field_entities.size() > 1) {
220 int dof_idx = 0;
221 for (auto dof : data.getFieldDofs()) {
222 auto &bit_ent = dof->getBitRefLevel();
223 if (!check(bitParentEnt, bitParentEntMask, bit_ent)) {
224 data.getIndices()[dof_idx] = -1;
225 data.getLocalIndices()[dof_idx] = -1;
226 data.getFieldData()[dof_idx] = 0;
227 }
228 ++dof_idx;
229 }
230 }
231
232#ifndef NDEBUG
233 if (verbosity >= NOISY) {
234 for (auto field_ent : field_entities) {
235 MOFEM_LOG("SELF", severityLevel)
236 << "Parent field entity bit: " << field_ent->getBitRefLevel()
237 << " " << *field_ent;
238 }
239 }
240#endif
241
242 if (field_entities.size()) {
243 boost::container::static_vector<EntityType, 128> ents_type;
244 ents_type.reserve(field_entities.size());
245 for (auto fe : field_entities)
246 ents_type.push_back(fe->getEntType());
247 std::sort(ents_type.begin(), ents_type.end());
248
249 auto end = std::unique(ents_type.begin(), ents_type.end());
250 for (auto it_t = ents_type.begin(); it_t != end; ++it_t)
251 CHKERR switch_of_dofs_children(
252 data, entities_field_data.dataOnEntities[*it_t]);
253 if (type == MBENTITYSET)
254 CHKERR switch_of_dofs_children(
255 data, entities_field_data.dataOnEntities[type]);
256 }
257
258 CHKERR set_child_data(data, child_data_meshset);
259 }
260 }
261
263 };
264
265 // iterate parents collect data
266 auto &bit_fe = getFEMethod()->numeredEntFiniteElementPtr->getBitRefLevel();
267 if (check(bitChild, bitChildMask, bit_fe)) {
268
269 if (verbosity >= VERBOSE) {
270 MOFEM_LOG("SELF", severityLevel) << "Child FE bit: " << bit_fe;
271 }
272
273 auto loop_parent_fe = [&]() {
275
276 parentElePtr->getOpPtrVector().back().doWorkRhsHook = do_work_parent_hook;
279
280#ifndef NDEBUG
281 auto &parent_gauss_pts = parentElePtr->gaussPts;
282 if (getGaussPts().size1() != parent_gauss_pts.size1()) {
283 MOFEM_LOG("SELF", Sev::error) << getGaussPts();
284 MOFEM_LOG("SELF", Sev::error) << parent_gauss_pts;
285 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
286 "Wrong number of weights");
287 }
288 if (getGaussPts().size2() != parent_gauss_pts.size2()) {
289 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
290 "Wrong number of integration points");
291 }
292#endif
293
295 };
296
297 CHKERR loop_parent_fe();
298 }
299
300 if (verbosity >= VERBOSE) {
301 MOFEM_LOG("SELF", severityLevel)
302 << "Number of meshset entities "
303 << entities_field_data.dataOnEntities[MBENTITYSET].size();
304 }
305
306 auto &data_on_meshset = entities_field_data.dataOnEntities[MBENTITYSET];
307 auto it = data_on_meshset.begin();
308
309 if (count_meshset_sides < data_on_meshset.size()) {
310 for (auto s = 0; s != count_meshset_sides; ++s)
311 ++it;
312 poolEntsVector.transfer(poolEntsVector.end(), it, data_on_meshset.end(),
313 data_on_meshset);
314 }
315
316 auto set_up_derivative_ent_data = [&](auto &entities_field_data,
317 auto &derivative_entities_field_data) {
319
321 using DerivedEntData = DerivedEntitiesFieldData::DerivedEntData;
322
324 auto &ents_data = entities_field_data.dataOnEntities[MBENTITYSET];
325 auto &dev_ents_data =
326 derivative_entities_field_data.dataOnEntities[MBENTITYSET];
327 dev_ents_data.clear();
328 for (auto c = 0; c < ents_data.size(); ++c) {
329 boost::shared_ptr<EntData> ent_data_ptr(data_ptr, &ents_data[c]);
330 dev_ents_data.push_back(new DerivedEntData(ent_data_ptr));
331 }
333 };
334
335 if (opParentType == OPSPACE) {
336 CHKERR set_up_derivative_ent_data(
337 entities_field_data,
338 *(getPtrFE()->getDerivedDataOnElementBySpaceArray()[approxSpace]));
339 }
340
341#ifndef NDEBUG
342 auto &side_table =
343 getPtrFE()->numeredEntFiniteElementPtr->getSideNumberTable();
344 for (auto &data : entities_field_data.dataOnEntities[MBENTITYSET]) {
345 for (auto dof : data.getFieldDofs()) {
346 auto &bit_dof = dof->getBitRefLevel();
347 if (check(bitParentEnt, bitParentEntMask, bit_dof)) {
348 auto ent = dof->getEnt();
349 if (side_table.find(ent) == side_table.end()) {
350 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
351 "Adjacency not found");
352 }
353 }
354 }
355 }
356#endif
357
359}
360
361} // namespace MoFEM
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:345
@ QUIET
Definition: definitions.h:208
@ NOISY
Definition: definitions.h:211
@ VERBOSE
Definition: definitions.h:209
FieldApproximationBase
approximation base
Definition: definitions.h:58
@ LASTBASE
Definition: definitions.h:69
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
FieldSpace
approximation spaces
Definition: definitions.h:82
@ NOSPACE
Definition: definitions.h:83
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
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
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 MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
FTensor::Index< 'm', SPACE_DIM > m
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
SeverityLevel
Severity levels.
Definition: LogManager.hpp:33
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:277
auto bit
set bit
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
constexpr auto field_name
base operator to do operations at Gauss Pt. level
Derived ata on single entity (This is passed as argument to DataOperator::doWork)
Data on single entity (This is passed as argument to DataOperator::doWork)
data structure for finite element entity
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
MoFEMErrorCode loopThis(const string &fe_name, ForcesAndSourcesCore *parent_fe, const int verb=QUIET, const LogManager::SeverityLevel sev=Sev::noisy)
User call this function to loop over parent elements. This function calls finite element with is oper...
MoFEMErrorCode loopParent(const string &fe_name, ForcesAndSourcesCore *parent_fe, const int verb=QUIET, const LogManager::SeverityLevel sev=Sev::noisy)
User call this function to loop over parent elements. This function calls finite element with is oper...
std::string getFEName() const
Get name of the element.
OpType
Controls loop over entities on element.
@ 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
structure to get information form mofem into EntitiesFieldData
auto & getDataOnElementBySpaceArray()
Get data on entities and space.
MoFEMErrorCode opRhs(EntitiesFieldData &data, const bool error_if_no_base=false)
boost::ptr_deque< EntitiesFieldData::EntData > poolEntsVector
boost::shared_ptr< ForcesAndSourcesCore > parentElePtr
OpAddParentEntData(std::string field_name, OpType op_parent_type, boost::shared_ptr< ForcesAndSourcesCore > parent_ele_ptr, BitRefLevel bit_child, BitRefLevel bit_child_mask, BitRefLevel bit_parent_ent, BitRefLevel bit_parent_ent_mask, int verb=QUIET, Sev sev=Sev::noisy)
Construct a new Op Add Parent Ent Data object.
boost::shared_ptr< ForcesAndSourcesCore > parentElePtr
boost::shared_ptr< ForcesAndSourcesCore > thisElePtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpRunParent(boost::shared_ptr< ForcesAndSourcesCore > parent_ele_ptr, BitRefLevel bit_parent, BitRefLevel bit_parent_mask, boost::shared_ptr< ForcesAndSourcesCore > this_ele_ptr, BitRefLevel bit_this, BitRefLevel bit_this_mask, int verb=QUIET, Sev sev=Sev::noisy)
Construct a new Op Run Parent object.