v0.13.0
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 /* MoFEM is free software: you can redistribute it and/or modify it under
7  * the terms of the GNU Lesser General Public License as published by the
8  * Free Software Foundation, either version 3 of the License, or (at your
9  * option) any later version.
10  *
11  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
12  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
14  * License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>
18  */
19 
20 namespace MoFEM {
21 
23 FieldBlas::query_interface(boost::typeindex::type_index type_index,
24  UnknownInterface **iface) const {
26  *iface = const_cast<FieldBlas *>(this);
28 }
29 
31  : cOre(const_cast<MoFEM::Core &>(core)), dEbug(false) {}
32 
35  const std::string field_name, Range *ents_ptr) {
36  const MoFEM::Interface &m_field = cOre;
37  auto fields_ptr = m_field.get_fields();
39 
40  auto x_fit = fields_ptr->get<FieldName_mi_tag>().find(field_name);
41  if (x_fit == fields_ptr->get<FieldName_mi_tag>().end()) {
42  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
43  "field < %s > not found, (top tip: check spelling)",
44  field_name.c_str());
45  }
46 
47  auto wrap_lambda_on_enties =
48  [lambda](boost::shared_ptr<FieldEntity> field_entity_ptr) {
50 
51  auto field_data = field_entity_ptr->getEntFieldData();
52  for (auto &v : field_data)
53  v = lambda(v);
54 
56  };
57 
58  CHKERR fieldLambdaOnEntities(wrap_lambda_on_enties, field_name, ents_ptr);
60 };
61 
64  const std::string field_name,
65  Range *ents_ptr) {
66  const MoFEM::Interface &m_field = cOre;
67  auto fields_ptr = m_field.get_fields();
68  auto field_ents = m_field.get_field_ents();
69  auto dofs_ptr = m_field.get_dofs();
71 
72  auto fit = fields_ptr->get<FieldName_mi_tag>().find(field_name);
73  if (fit == fields_ptr->get<FieldName_mi_tag>().end()) {
74  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
75  "field < %s > not found, (top tip: check spelling)",
76  field_name.c_str());
77  }
78 
79  // End of iterator for y_field
80  const auto bit_number = (*fit)->getBitNumber();
81 
82  auto loop_ents = [&](const EntityHandle f, const EntityHandle s) {
84 
85  const auto lo_uid = FieldEntity::getLoLocalEntityBitNumber(bit_number, f);
86  const auto hi_uid = FieldEntity::getHiLocalEntityBitNumber(bit_number, s);
87 
88  // Start of iterator for x_field
89  auto eit = field_ents->get<Unique_mi_tag>().lower_bound(lo_uid);
90  // End of iterator for x_field
91  auto eit_hi = field_ents->get<Unique_mi_tag>().upper_bound(hi_uid);
92 
93  for (; eit != eit_hi; ++eit) {
94  CHKERR lambda(*eit);
95  }
96 
98  };
99 
100  if (ents_ptr) {
101  for (auto p = ents_ptr->const_pair_begin(); p != ents_ptr->const_pair_end();
102  ++p) {
103  CHKERR loop_ents(p->first, p->second);
104  }
105  } else {
106  // we are looping from the very first possible entity handle (MBVERTEX) to
107  // the very last (MBENTITYSET)
108  CHKERR loop_ents(get_id_for_min_type<MBVERTEX>(),
109  get_id_for_max_type<MBENTITYSET>());
110  }
111 
113 }
114 
116  FieldBlas::TwoFieldFunctionOnValues lambda, const std::string field_name_x,
117  const std::string field_name_y, bool error_if_missing, Range *ents_ptr) {
118  const MoFEM::Interface &m_field = cOre;
119  auto fields_ptr = m_field.get_fields();
121 
122  auto x_fit = fields_ptr->get<FieldName_mi_tag>().find(field_name_x);
123  if (x_fit == fields_ptr->get<FieldName_mi_tag>().end()) {
124  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
125  "x field < %s > not found, (top tip: check spelling)",
126  field_name_x.c_str());
127  }
128  auto y_fit = fields_ptr->get<FieldName_mi_tag>().find(field_name_y);
129  if (y_fit == fields_ptr->get<FieldName_mi_tag>().end()) {
130  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
131  "y field < %s > not found, (top tip: check spelling)",
132  field_name_y.c_str());
133  }
134  if ((*x_fit)->getSpace() != (*y_fit)->getSpace()) {
135  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
136  "space for field < %s > and field <%s> are not compatible for "
137  "fieldblas",
138  field_name_x.c_str(), field_name_y.c_str());
139  }
140  if ((*x_fit)->getNbOfCoeffs() != (*y_fit)->getNbOfCoeffs()) {
141  SETERRQ2(
142  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
143  "rank for field < %s > and field <%s> are not compatible for fieldblas",
144  field_name_x.c_str(), field_name_y.c_str());
145  }
146 
147  auto wrap_lambda_on_enties =
148  [lambda](boost::shared_ptr<FieldEntity> y_field_entity_ptr,
149  const boost::shared_ptr<FieldEntity> x_field_entity_ptr) {
151 
152  auto x_field_data = x_field_entity_ptr->getEntFieldData();
153  auto y_field_data = y_field_entity_ptr->getEntFieldData();
154  const auto size_x = x_field_data.size();
155  const auto size_y = y_field_data.size();
156 
157  size_t dd = 0;
158  for (; dd != std::min(size_x, size_y); ++dd)
159  y_field_data[dd] = lambda(y_field_data[dd], x_field_data[dd]);
160  for (; dd < size_y; ++dd)
161  y_field_data[dd] = 0;
162 
164  };
165 
166  CHKERR fieldLambdaOnEntities(wrap_lambda_on_enties, field_name_x,
167  field_name_y, error_if_missing, ents_ptr);
169 };
170 
173  const std::string field_name_x,
174  const std::string field_name_y,
175  bool error_if_missing, Range *ents_ptr) {
176  const MoFEM::Interface &m_field = cOre;
177  auto fields_ptr = m_field.get_fields();
178  auto field_ents = m_field.get_field_ents();
179  auto dofs_ptr = m_field.get_dofs();
181 
182  auto x_fit = fields_ptr->get<FieldName_mi_tag>().find(field_name_x);
183  if (x_fit == fields_ptr->get<FieldName_mi_tag>().end()) {
184  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
185  "x field < %s > not found, (top tip: check spelling)",
186  field_name_x.c_str());
187  }
188  auto y_fit = fields_ptr->get<FieldName_mi_tag>().find(field_name_y);
189  if (y_fit == fields_ptr->get<FieldName_mi_tag>().end()) {
190  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
191  "y field < %s > not found, (top tip: check spelling)",
192  field_name_y.c_str());
193  }
194 
195  // End of iterator for y_field
196  const auto x_bit_number = (*x_fit)->getBitNumber();
197  const auto y_bit_number = (*y_fit)->getBitNumber();
198  const auto y_eit_hi = field_ents->get<Unique_mi_tag>().upper_bound(
199  FieldEntity::getHiBitNumberUId(y_bit_number));
200 
201  auto loop_ents = [&](const EntityHandle f, const EntityHandle s) {
203 
204  const auto x_lo_uid =
206  const auto x_hi_uid =
208 
209  // Start of iterator for x_field
210  auto x_eit = field_ents->get<Unique_mi_tag>().lower_bound(x_lo_uid);
211  // End of iterator for x_field
212  auto x_eit_hi = field_ents->get<Unique_mi_tag>().upper_bound(x_hi_uid);
213 
214  for (; x_eit != x_eit_hi;) {
215 
216  const auto y_lo_uid = FieldEntity::getLocalUniqueIdCalculate(
217  (*y_fit)->getBitNumber(), (*x_eit)->getEnt());
218  auto y_eit = field_ents->get<Unique_mi_tag>().find(y_lo_uid);
219 
220  if (y_eit == field_ents->end()) {
221 
222  if (error_if_missing) {
223  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
224  "Missing entity in y_field.");
225  }
226 
227  ++x_eit;
228 
229  } else {
230 
231  auto check = [&]() {
232  if (x_eit == x_eit_hi)
233  return false;
234  if (y_eit == y_eit_hi)
235  return false;
236  if ((*y_eit)->getEnt() != (*x_eit)->getEnt())
237  return false;
238  return true;
239  };
240 
241  do {
242 
243  CHKERR lambda(*y_eit, *x_eit);
244 
245  ++x_eit;
246  ++y_eit;
247 
248  } while (check());
249  }
250  }
251 
253  };
254 
255  if (ents_ptr) {
256  for (auto p = ents_ptr->const_pair_begin(); p != ents_ptr->const_pair_end();
257  ++p) {
258  CHKERR loop_ents(p->first, p->second);
259  }
260  } else {
261  // we are looping from the very first possible entity handle (MBVERTEX) to
262  // the very last (MBENTITYSET)
263  CHKERR loop_ents(get_id_for_min_type<MBVERTEX>(),
264  get_id_for_max_type<MBENTITYSET>());
265  }
266 
268 }
269 
271  const std::string field_name_x,
272  const std::string field_name_y,
273  bool error_if_missing,
274  bool create_if_missing) {
275  if (create_if_missing)
276  MOFEM_LOG("SELF", Sev::noisy)
277  << "Option create_if_missing is set to true but have no effect";
278  return fieldLambdaOnValues(lambda, field_name_x, field_name_y,
279  error_if_missing);
280 }
281 
283  const std::string field_name_x,
284  const std::string field_name_y,
285  bool error_if_missing, Range *ents_ptr) {
287  auto axpy = [alpha](const double fy, const double fx) {
288  return fy + alpha * fx;
289  };
290  CHKERR fieldLambdaOnValues(axpy, field_name_x, field_name_y, error_if_missing,
291  ents_ptr);
293 }
294 
296  const std::string field_name_x,
297  const std::string field_name_y,
298  bool error_if_missing,
299  bool create_if_missing) {
300  return fieldAxpy(alpha, field_name_x, field_name_y, error_if_missing);
301 }
302 
304  const std::string field_name_x,
305  const std::string field_name_y,
306  bool error_if_missing, Range *ents_ptr) {
308  auto copy = [alpha](const double fy, const double fx) { return alpha * fx; };
309  CHKERR fieldLambdaOnValues(copy, field_name_x, field_name_y, error_if_missing,
310  ents_ptr);
312 }
313 
315  const std::string field_name_x,
316  const std::string field_name_y,
317  bool error_if_missing,
318  bool create_if_missing) {
319  return fieldCopy(alpha, field_name_x, field_name_y, error_if_missing);
320 }
321 
323  const std::string field_name,
324  Range *ents_ptr) {
326 
327  auto scale = [alpha](const double v) { return v * alpha; };
328  CHKERR fieldLambdaOnValues(scale, field_name, ents_ptr);
329 
331 }
332 
334  const std::string field_name,
335  Range *sub_verts) {
336  const MoFEM::Interface &m_field = cOre;
338 
339  EntityHandle meshset = m_field.get_field_meshset(field_name);
340  Range verts;
341  CHKERR m_field.get_moab().get_entities_by_type(meshset, MBVERTEX, verts,
342  true);
343  if (sub_verts)
344  verts = intersect(*sub_verts, verts);
345 
346  struct LambdaMethod : EntityMethod {
347  LambdaMethod(MoFEM::Interface &m_field, Range &verts,
349  : EntityMethod(), mField(m_field), verts(verts), lambda(lambda),
350  count(0), total(0) {}
351  MoFEMErrorCode preProcess() {
352  vit = verts.begin();
353  return 0;
354  }
355  MoFEMErrorCode operator()() {
357  if (*vit != entPtr->getEnt())
358  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
359  "Inconsistent entity %ld != %ld", *vit, entPtr->getEnt());
360  if (!count)
361  CHKERR mField.get_moab().coords_iterate(vit, verts.end(), xCoord,
362  yCoord, zCoord, count);
363  CHKERR lambda(entPtr->getEntFieldData(), xCoord, yCoord, zCoord);
364  ++xCoord;
365  ++yCoord;
366  ++zCoord;
367  ++vit;
368  ++total;
369  --count;
371  }
372  MoFEMErrorCode postProcess() {
374  if (total != verts.size())
375  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
376  "Inconsistent number of vertices in field meshset and in the "
377  "field %d != %d",
378  total, verts.size());
380  }
381 
382  private:
383  MoFEM::Interface &mField;
384  Range &verts;
386  int count;
387  int total;
388  Range::iterator vit;
389  double *xCoord;
390  double *yCoord;
391  double *zCoord;
392  };
393 
394  LambdaMethod lambda_method(const_cast<MoFEM::Interface &>(m_field), verts,
395  lambda);
396  CHKERR const_cast<MoFEM::Interface &>(m_field).loop_entities(
397  field_name, lambda_method, &verts, QUIET);
399 }
400 
401 MoFEMErrorCode FieldBlas::setField(const double val, const EntityType type,
402  const std::string field_name) {
403  const MoFEM::Interface &m_field = cOre;
404  auto dofs_ptr = m_field.get_dofs();
406 
407  const auto bit_number = m_field.get_field_bit_number(field_name);
408  const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
409  bit_number, get_id_for_min_type(type));
410  const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
411  bit_number, get_id_for_max_type(type));
412 
413  auto dit = dofs_ptr->get<Unique_mi_tag>().lower_bound(lo_uid);
414  auto hi_dit = dofs_ptr->get<Unique_mi_tag>().upper_bound(hi_uid);
415 
416  for (; dit != hi_dit; dit++)
417  (*dit)->getFieldData() = val;
418 
420 }
421 
422 MoFEMErrorCode FieldBlas::setField(const double val, const EntityType type,
423  const Range &ents,
424  const std::string field_name) {
425  const MoFEM::Interface &m_field = cOre;
426  auto dofs_ptr = m_field.get_dofs();
428 
429  const auto bit_number = m_field.get_field_bit_number(field_name);
430  const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
431  bit_number, get_id_for_min_type(type));
432  const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
433  bit_number, get_id_for_max_type(type));
434 
435  auto dit = dofs_ptr->get<Unique_mi_tag>().lower_bound(lo_uid);
436  auto hi_dit = dofs_ptr->get<Unique_mi_tag>().upper_bound(hi_uid);
437 
438  EntityHandle ent, last = 0;
439  bool cont = true;
440  for (; dit != hi_dit; dit++) {
441  ent = (*dit)->getEnt();
442  if (ent != last) {
443  if (ents.find(ent) == ents.end()) {
444  cont = true;
445  } else {
446  cont = false;
447  }
448  last = ent;
449  }
450  if (!cont)
451  (*dit)->getFieldData() = val;
452  }
454 }
455 
457  const std::string field_name) {
458  const MoFEM::Interface &m_field = cOre;
459  auto fields_ptr = m_field.get_fields();
460  auto dofs_ptr = m_field.get_dofs();
462  auto fit = fields_ptr->get<FieldName_mi_tag>().find(field_name);
463  if (fit == fields_ptr->get<FieldName_mi_tag>().end()) {
464  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
465  " field < %s > not found, (top tip: check spelling)",
466  field_name.c_str());
467  }
468 
469  auto dit = dofs_ptr->get<Unique_mi_tag>().lower_bound(
470  FieldEntity::getLoBitNumberUId((*fit)->getBitNumber()));
471  auto hi_dit = dofs_ptr->get<Unique_mi_tag>().upper_bound(
472  FieldEntity::getHiBitNumberUId((*fit)->getBitNumber()));
473  for (; dit != hi_dit; dit++) {
474  (*dit)->getFieldData() = val;
475  }
477 }
478 
479 } // namespace MoFEM
static Index< 'p', 3 > p
constexpr double alpha
@ QUIET
Definition: definitions.h:221
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
constexpr double lambda
virtual const FieldEntity_multiIndex * get_field_ents() const =0
Get the field ents object.
virtual const Field_multiIndex * get_fields() const =0
Get the fields object.
virtual const DofEntity_multiIndex * get_dofs() const =0
Get the dofs object.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:311
double v
phase velocity of light in medium (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)
Definition: ddTensor0.hpp:33
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
EntityHandle get_id_for_max_type()
EntityHandle get_id_for_min_type()
double scale
Definition: plastic.cpp:103
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
virtual EntityHandle get_field_meshset(const std::string name) const =0
get field meshset
virtual moab::Interface & get_moab()=0
Core (interface) class.
Definition: Core.hpp:92
Deprecated interface functions.
Data structure to exchange data between mofem and User Loop Methods on entities.
Basic algebra on fields.
Definition: FieldBlas.hpp:31
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:270
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:282
boost::function< MoFEMErrorCode(boost::shared_ptr< FieldEntity > field_entity_ptr)> OneFieldFunctionOnEntities
Definition: FieldBlas.hpp:55
MoFEMErrorCode fieldLambdaOnEntities(OneFieldFunctionOnEntities lambda, const std::string field_name, Range *ents_ptr=nullptr)
filed lambda
Definition: FieldBlas.cpp:63
boost::function< double(const double y_field_value_reference, const double x_field_value)> TwoFieldFunctionOnValues
Definition: FieldBlas.hpp:63
MoFEMErrorCode setVertexDofs(VertexCoordsFunction lambda, const std::string field_name, Range *verts=nullptr)
Set DOFs on vertices using user function.
Definition: FieldBlas.cpp:333
boost::function< MoFEMErrorCode(VectorAdaptor &&field_data, double *xcoord, double *ycoord, double *zcoord)> VertexCoordsFunction
Definition: FieldBlas.hpp:292
MoFEMErrorCode setField(const double val, const EntityType type, const std::string field_name)
Definition: FieldBlas.cpp:401
boost::function< double(const double field_val)> OneFieldFunctionOnValues
function to set a field value
Definition: FieldBlas.hpp:48
MoFEMErrorCode fieldScale(const double alpha, const std::string field_name, Range *ents_ptr=nullptr)
field scale
Definition: FieldBlas.cpp:322
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Definition: FieldBlas.cpp:23
MoFEMErrorCode fieldLambdaOnValues(OneFieldFunctionOnValues lambda, const std::string field_name, Range *ents_ptr=nullptr)
filed lambda
Definition: FieldBlas.cpp:34
FieldBlas(const MoFEM::Core &core)
Definition: FieldBlas.cpp:30
const MoFEM::Interface & cOre
Definition: FieldBlas.hpp:36
boost::function< MoFEMErrorCode(boost::shared_ptr< FieldEntity > y_field_entity_ptr, const boost::shared_ptr< FieldEntity > x_field_entity_ptr)> TwoFieldFunctionOnEntities
Definition: FieldBlas.hpp:72
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:303
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