v0.13.1
Loading...
Searching...
No Matches
FieldBlas.cpp
Go to the documentation of this file.
1/** \file FieldBlas.cpp
2 * \brief Managing complexities for problem
3 * \ingroup mofem_section_manager
4 */
5
6
7namespace MoFEM {
8
10FieldBlas::query_interface(boost::typeindex::type_index type_index,
11 UnknownInterface **iface) const {
13 *iface = const_cast<FieldBlas *>(this);
15}
16
18 : cOre(const_cast<MoFEM::Core &>(core)), dEbug(false) {}
19
22 const std::string field_name, Range *ents_ptr) {
23 const MoFEM::Interface &m_field = cOre;
24 auto fields_ptr = m_field.get_fields();
26
27 auto x_fit = fields_ptr->get<FieldName_mi_tag>().find(field_name);
28 if (x_fit == fields_ptr->get<FieldName_mi_tag>().end()) {
29 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
30 "field < %s > not found, (top tip: check spelling)",
31 field_name.c_str());
32 }
33
34 auto wrap_lambda_on_enties =
35 [lambda](boost::shared_ptr<FieldEntity> field_entity_ptr) {
37
38 auto field_data = field_entity_ptr->getEntFieldData();
39 for (auto &v : field_data)
40 v = lambda(v);
41
43 };
44
45 CHKERR fieldLambdaOnEntities(wrap_lambda_on_enties, field_name, ents_ptr);
47};
48
51 const std::string field_name,
52 Range *ents_ptr) {
53 const MoFEM::Interface &m_field = cOre;
54 auto fields_ptr = m_field.get_fields();
55 auto field_ents = m_field.get_field_ents();
56 auto dofs_ptr = m_field.get_dofs();
58
59 auto fit = fields_ptr->get<FieldName_mi_tag>().find(field_name);
60 if (fit == fields_ptr->get<FieldName_mi_tag>().end()) {
61 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
62 "field < %s > not found, (top tip: check spelling)",
63 field_name.c_str());
64 }
65
66 // End of iterator for y_field
67 const auto bit_number = (*fit)->getBitNumber();
68
69 auto loop_ents = [&](const EntityHandle f, const EntityHandle s) {
71
72 const auto lo_uid = FieldEntity::getLoLocalEntityBitNumber(bit_number, f);
73 const auto hi_uid = FieldEntity::getHiLocalEntityBitNumber(bit_number, s);
74
75 // Start of iterator for x_field
76 auto eit = field_ents->get<Unique_mi_tag>().lower_bound(lo_uid);
77 // End of iterator for x_field
78 auto eit_hi = field_ents->get<Unique_mi_tag>().upper_bound(hi_uid);
79
80 for (; eit != eit_hi; ++eit) {
81 CHKERR lambda(*eit);
82 }
83
85 };
86
87 if (ents_ptr) {
88 for (auto p = ents_ptr->const_pair_begin(); p != ents_ptr->const_pair_end();
89 ++p) {
90 CHKERR loop_ents(p->first, p->second);
91 }
92 } else {
93 // we are looping from the very first possible entity handle (MBVERTEX) to
94 // the very last (MBENTITYSET)
95 CHKERR loop_ents(get_id_for_min_type<MBVERTEX>(),
96 get_id_for_max_type<MBENTITYSET>());
97 }
98
100}
101
103 FieldBlas::TwoFieldFunctionOnValues lambda, const std::string field_name_x,
104 const std::string field_name_y, bool error_if_missing, Range *ents_ptr) {
105 const MoFEM::Interface &m_field = cOre;
106 auto fields_ptr = m_field.get_fields();
108
109 auto x_fit = fields_ptr->get<FieldName_mi_tag>().find(field_name_x);
110 if (x_fit == fields_ptr->get<FieldName_mi_tag>().end()) {
111 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
112 "x field < %s > not found, (top tip: check spelling)",
113 field_name_x.c_str());
114 }
115 auto y_fit = fields_ptr->get<FieldName_mi_tag>().find(field_name_y);
116 if (y_fit == fields_ptr->get<FieldName_mi_tag>().end()) {
117 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
118 "y field < %s > not found, (top tip: check spelling)",
119 field_name_y.c_str());
120 }
121 if ((*x_fit)->getSpace() != (*y_fit)->getSpace()) {
122 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
123 "space for field < %s > and field <%s> are not compatible for "
124 "fieldblas",
125 field_name_x.c_str(), field_name_y.c_str());
126 }
127 if ((*x_fit)->getNbOfCoeffs() != (*y_fit)->getNbOfCoeffs()) {
128 SETERRQ2(
129 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
130 "rank for field < %s > and field <%s> are not compatible for fieldblas",
131 field_name_x.c_str(), field_name_y.c_str());
132 }
133
134 auto wrap_lambda_on_enties =
135 [lambda](boost::shared_ptr<FieldEntity> y_field_entity_ptr,
136 const boost::shared_ptr<FieldEntity> x_field_entity_ptr) {
138
139 auto x_field_data = x_field_entity_ptr->getEntFieldData();
140 auto y_field_data = y_field_entity_ptr->getEntFieldData();
141 const auto size_x = x_field_data.size();
142 const auto size_y = y_field_data.size();
143
144 size_t dd = 0;
145 for (; dd != std::min(size_x, size_y); ++dd)
146 y_field_data[dd] = lambda(y_field_data[dd], x_field_data[dd]);
147 for (; dd < size_y; ++dd)
148 y_field_data[dd] = 0;
149
151 };
152
153 CHKERR fieldLambdaOnEntities(wrap_lambda_on_enties, field_name_x,
154 field_name_y, error_if_missing, ents_ptr);
156};
157
160 const std::string field_name_x,
161 const std::string field_name_y,
162 bool error_if_missing, Range *ents_ptr) {
163 const MoFEM::Interface &m_field = cOre;
164 auto fields_ptr = m_field.get_fields();
165 auto field_ents = m_field.get_field_ents();
166 auto dofs_ptr = m_field.get_dofs();
168
169 auto x_fit = fields_ptr->get<FieldName_mi_tag>().find(field_name_x);
170 if (x_fit == fields_ptr->get<FieldName_mi_tag>().end()) {
171 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
172 "x field < %s > not found, (top tip: check spelling)",
173 field_name_x.c_str());
174 }
175 auto y_fit = fields_ptr->get<FieldName_mi_tag>().find(field_name_y);
176 if (y_fit == fields_ptr->get<FieldName_mi_tag>().end()) {
177 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
178 "y field < %s > not found, (top tip: check spelling)",
179 field_name_y.c_str());
180 }
181
182 // End of iterator for y_field
183 const auto x_bit_number = (*x_fit)->getBitNumber();
184 const auto y_bit_number = (*y_fit)->getBitNumber();
185 const auto y_eit_hi = field_ents->get<Unique_mi_tag>().upper_bound(
186 FieldEntity::getHiBitNumberUId(y_bit_number));
187
188 auto loop_ents = [&](const EntityHandle f, const EntityHandle s) {
190
191 const auto x_lo_uid =
193 const auto x_hi_uid =
195
196 // Start of iterator for x_field
197 auto x_eit = field_ents->get<Unique_mi_tag>().lower_bound(x_lo_uid);
198 // End of iterator for x_field
199 auto x_eit_hi = field_ents->get<Unique_mi_tag>().upper_bound(x_hi_uid);
200
201 for (; x_eit != x_eit_hi;) {
202
203 const auto y_lo_uid = FieldEntity::getLocalUniqueIdCalculate(
204 (*y_fit)->getBitNumber(), (*x_eit)->getEnt());
205 auto y_eit = field_ents->get<Unique_mi_tag>().find(y_lo_uid);
206
207 if (y_eit == field_ents->end()) {
208
209 if (error_if_missing) {
210 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
211 "Missing entity in y_field.");
212 }
213
214 ++x_eit;
215
216 } else {
217
218 auto check = [&]() {
219 if (x_eit == x_eit_hi)
220 return false;
221 if (y_eit == y_eit_hi)
222 return false;
223 if ((*y_eit)->getEnt() != (*x_eit)->getEnt())
224 return false;
225 return true;
226 };
227
228 do {
229
230 CHKERR lambda(*y_eit, *x_eit);
231
232 ++x_eit;
233 ++y_eit;
234
235 } while (check());
236 }
237 }
238
240 };
241
242 if (ents_ptr) {
243 for (auto p = ents_ptr->const_pair_begin(); p != ents_ptr->const_pair_end();
244 ++p) {
245 CHKERR loop_ents(p->first, p->second);
246 }
247 } else {
248 // we are looping from the very first possible entity handle (MBVERTEX) to
249 // the very last (MBENTITYSET)
250 CHKERR loop_ents(get_id_for_min_type<MBVERTEX>(),
251 get_id_for_max_type<MBENTITYSET>());
252 }
253
255}
256
258 const std::string field_name_x,
259 const std::string field_name_y,
260 bool error_if_missing,
261 bool create_if_missing) {
262 if (create_if_missing)
263 MOFEM_LOG("SELF", Sev::noisy)
264 << "Option create_if_missing is set to true but have no effect";
265 return fieldLambdaOnValues(lambda, field_name_x, field_name_y,
266 error_if_missing);
267}
268
270 const std::string field_name_x,
271 const std::string field_name_y,
272 bool error_if_missing, Range *ents_ptr) {
274 auto axpy = [alpha](const double fy, const double fx) {
275 return fy + alpha * fx;
276 };
277 CHKERR fieldLambdaOnValues(axpy, field_name_x, field_name_y, error_if_missing,
278 ents_ptr);
280}
281
283 const std::string field_name_x,
284 const std::string field_name_y,
285 bool error_if_missing,
286 bool create_if_missing) {
287 return fieldAxpy(alpha, field_name_x, field_name_y, error_if_missing);
288}
289
291 const std::string field_name_x,
292 const std::string field_name_y,
293 bool error_if_missing, Range *ents_ptr) {
295 auto copy = [alpha](const double fy, const double fx) { return alpha * fx; };
296 CHKERR fieldLambdaOnValues(copy, field_name_x, field_name_y, error_if_missing,
297 ents_ptr);
299}
300
302 const std::string field_name_x,
303 const std::string field_name_y,
304 bool error_if_missing,
305 bool create_if_missing) {
306 return fieldCopy(alpha, field_name_x, field_name_y, error_if_missing);
307}
308
310 const std::string field_name,
311 Range *ents_ptr) {
313
314 auto scale = [alpha](const double v) { return v * alpha; };
316
318}
319
321 const std::string field_name,
322 Range *sub_verts) {
323 const MoFEM::Interface &m_field = cOre;
325
326 EntityHandle meshset = m_field.get_field_meshset(field_name);
327 Range verts;
328 CHKERR m_field.get_moab().get_entities_by_type(meshset, MBVERTEX, verts,
329 true);
330 if (sub_verts)
331 verts = intersect(*sub_verts, verts);
332
333 struct LambdaMethod : EntityMethod {
334 LambdaMethod(MoFEM::Interface &m_field, Range &verts,
336 : EntityMethod(), mField(m_field), verts(verts), lambda(lambda),
337 count(0), total(0) {}
338 MoFEMErrorCode preProcess() {
339 vit = verts.begin();
340 return 0;
341 }
342 MoFEMErrorCode operator()() {
344 if (*vit != entPtr->getEnt())
345 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
346 "Inconsistent entity %ld != %ld", *vit, entPtr->getEnt());
347 if (!count)
348 CHKERR mField.get_moab().coords_iterate(vit, verts.end(), xCoord,
349 yCoord, zCoord, count);
350 CHKERR lambda(entPtr->getEntFieldData(), xCoord, yCoord, zCoord);
351 ++xCoord;
352 ++yCoord;
353 ++zCoord;
354 ++vit;
355 ++total;
356 --count;
358 }
359 MoFEMErrorCode postProcess() {
361 if (total != verts.size())
362 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
363 "Inconsistent number of vertices in field meshset and in the "
364 "field %d != %d",
365 total, verts.size());
367 }
368
369 private:
370 MoFEM::Interface &mField;
371 Range &verts;
373 int count;
374 int total;
375 Range::iterator vit;
376 double *xCoord;
377 double *yCoord;
378 double *zCoord;
379 };
380
381 LambdaMethod lambda_method(const_cast<MoFEM::Interface &>(m_field), verts,
382 lambda);
383 CHKERR const_cast<MoFEM::Interface &>(m_field).loop_entities(
384 field_name, lambda_method, &verts, QUIET);
386}
387
388MoFEMErrorCode FieldBlas::setField(const double val, const EntityType type,
389 const std::string field_name) {
390 const MoFEM::Interface &m_field = cOre;
391 auto dofs_ptr = m_field.get_dofs();
393
394 const auto bit_number = m_field.get_field_bit_number(field_name);
395 const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
396 bit_number, get_id_for_min_type(type));
397 const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
398 bit_number, get_id_for_max_type(type));
399
400 auto dit = dofs_ptr->get<Unique_mi_tag>().lower_bound(lo_uid);
401 auto hi_dit = dofs_ptr->get<Unique_mi_tag>().upper_bound(hi_uid);
402
403 for (; dit != hi_dit; dit++)
404 (*dit)->getFieldData() = val;
405
407}
408
409MoFEMErrorCode FieldBlas::setField(const double val, const EntityType type,
410 const Range &ents,
411 const std::string field_name) {
412 const MoFEM::Interface &m_field = cOre;
413 auto dofs_ptr = m_field.get_dofs();
415
416 const auto bit_number = m_field.get_field_bit_number(field_name);
417 const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
418 bit_number, get_id_for_min_type(type));
419 const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
420 bit_number, get_id_for_max_type(type));
421
422 auto dit = dofs_ptr->get<Unique_mi_tag>().lower_bound(lo_uid);
423 auto hi_dit = dofs_ptr->get<Unique_mi_tag>().upper_bound(hi_uid);
424
425 EntityHandle ent, last = 0;
426 bool cont = true;
427 for (; dit != hi_dit; dit++) {
428 ent = (*dit)->getEnt();
429 if (ent != last) {
430 if (ents.find(ent) == ents.end()) {
431 cont = true;
432 } else {
433 cont = false;
434 }
435 last = ent;
436 }
437 if (!cont)
438 (*dit)->getFieldData() = val;
439 }
441}
442
444 const std::string field_name) {
445 const MoFEM::Interface &m_field = cOre;
446 auto fields_ptr = m_field.get_fields();
447 auto dofs_ptr = m_field.get_dofs();
449 auto fit = fields_ptr->get<FieldName_mi_tag>().find(field_name);
450 if (fit == fields_ptr->get<FieldName_mi_tag>().end()) {
451 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
452 " field < %s > not found, (top tip: check spelling)",
453 field_name.c_str());
454 }
455
456 auto dit = dofs_ptr->get<Unique_mi_tag>().lower_bound(
457 FieldEntity::getLoBitNumberUId((*fit)->getBitNumber()));
458 auto hi_dit = dofs_ptr->get<Unique_mi_tag>().upper_bound(
459 FieldEntity::getHiBitNumberUId((*fit)->getBitNumber()));
460 for (; dit != hi_dit; dit++) {
461 (*dit)->getFieldData() = val;
462 }
464}
465
466} // namespace MoFEM
static Index< 'p', 3 > p
@ QUIET
Definition: definitions.h:208
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#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
#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
constexpr double lambda
virtual const FieldEntity_multiIndex * get_field_ents() const =0
Get the field ents object.
virtual const DofEntity_multiIndex * get_dofs() const =0
Get the dofs object.
virtual const Field_multiIndex * get_fields() const =0
Get the fields object.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
const double v
phase velocity of light in medium (cm/ns)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
EntityHandle get_id_for_max_type()
EntityHandle get_id_for_min_type()
double scale
Definition: plastic.cpp:94
constexpr auto field_name
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
virtual moab::Interface & get_moab()=0
virtual EntityHandle get_field_meshset(const std::string name) const =0
get field meshset
Core (interface) class.
Definition: Core.hpp:82
Deprecated interface functions.
Data structure to exchange data between mofem and User Loop Methods on entities.
Basic algebra on fields.
Definition: FieldBlas.hpp:21
DEPRECATED MoFEMErrorCode fieldLambda(TwoFieldFunctionOnValues lambda, const std::string field_name_x, const std::string field_name_y, bool error_if_missing, bool create_if_missing)
Definition: FieldBlas.cpp:257
MoFEMErrorCode fieldAxpy(const double alpha, const std::string field_name_x, const std::string field_name_y, bool error_if_missing=false, Range *ents_ptr=nullptr)
axpy fields
Definition: FieldBlas.cpp:269
boost::function< MoFEMErrorCode(boost::shared_ptr< FieldEntity > field_entity_ptr)> OneFieldFunctionOnEntities
Definition: FieldBlas.hpp:45
MoFEMErrorCode fieldLambdaOnEntities(OneFieldFunctionOnEntities lambda, const std::string field_name, Range *ents_ptr=nullptr)
filed lambda
Definition: FieldBlas.cpp:50
boost::function< double(const double y_field_value_reference, const double x_field_value)> TwoFieldFunctionOnValues
Definition: FieldBlas.hpp:53
MoFEMErrorCode setVertexDofs(VertexCoordsFunction lambda, const std::string field_name, Range *verts=nullptr)
Set DOFs on vertices using user function.
Definition: FieldBlas.cpp:320
boost::function< MoFEMErrorCode(VectorAdaptor &&field_data, double *xcoord, double *ycoord, double *zcoord)> VertexCoordsFunction
Definition: FieldBlas.hpp:282
MoFEMErrorCode setField(const double val, const EntityType type, const std::string field_name)
Definition: FieldBlas.cpp:388
boost::function< double(const double field_val)> OneFieldFunctionOnValues
function to set a field value
Definition: FieldBlas.hpp:38
MoFEMErrorCode fieldScale(const double alpha, const std::string field_name, Range *ents_ptr=nullptr)
field scale
Definition: FieldBlas.cpp:309
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Definition: FieldBlas.cpp:10
MoFEMErrorCode fieldLambdaOnValues(OneFieldFunctionOnValues lambda, const std::string field_name, Range *ents_ptr=nullptr)
filed lambda
Definition: FieldBlas.cpp:21
FieldBlas(const MoFEM::Core &core)
Definition: FieldBlas.cpp:17
const MoFEM::Interface & cOre
Definition: FieldBlas.hpp:26
boost::function< MoFEMErrorCode(boost::shared_ptr< FieldEntity > y_field_entity_ptr, const boost::shared_ptr< FieldEntity > x_field_entity_ptr)> TwoFieldFunctionOnEntities
Definition: FieldBlas.hpp:62
MoFEMErrorCode fieldCopy(const double alpha, const std::string field_name_x, const std::string field_name_y, bool error_if_missing=false, Range *ents_ptr=nullptr)
copy and scale fields
Definition: FieldBlas.cpp:290
static UId getHiBitNumberUId(const FieldBitNumber bit_number)
UId getLocalUniqueIdCalculate()
Get the Local Unique Id Calculate object.
static UId getLoBitNumberUId(const FieldBitNumber bit_number)
static UId getLoLocalEntityBitNumber(const char bit_number, const EntityHandle ent)
static UId getHiLocalEntityBitNumber(const char bit_number, const EntityHandle ent)
MultiIndex Tag for field name.
base class for all interface classes