v0.9.2
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  UnknownInterface **iface) const {
25  *iface = NULL;
26  if (uuid == IDD_MOFEMFieldBlas) {
27  *iface = const_cast<FieldBlas *>(this);
29  }
30  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "unknown interface");
32 }
33 
35  : cOre(const_cast<MoFEM::Core &>(core)), dEbug(false) {}
37 
39  const std::string field_name_x,
40  const std::string field_name_y,
41  bool error_if_missing,
42  bool creat_if_missing) {
43  const MoFEM::Interface &m_field = cOre;
44  auto fields_ptr = m_field.get_fields();
45  auto field_ents = m_field.get_field_ents();
46  auto dofs_ptr = m_field.get_dofs();
48 
49  auto x_fit = fields_ptr->get<FieldName_mi_tag>().find(field_name_x);
50  if (x_fit == fields_ptr->get<FieldName_mi_tag>().end()) {
51  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
52  "x field < %s > not found, (top tip: check spelling)",
53  field_name_x.c_str());
54  }
55  auto y_fit = fields_ptr->get<FieldName_mi_tag>().find(field_name_y);
56  if (y_fit == fields_ptr->get<FieldName_mi_tag>().end()) {
57  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
58  "y field < %s > not found, (top tip: check spelling)",
59  field_name_y.c_str());
60  }
61  if ((*x_fit)->getSpace() != (*y_fit)->getSpace()) {
62  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
63  "space for field < %s > and field <%s> are not compatible",
64  field_name_x.c_str(), field_name_y.c_str());
65  }
66  if ((*x_fit)->getNbOfCoeffs() != (*y_fit)->getNbOfCoeffs()) {
67  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
68  "rank for field < %s > and field <%s> are not compatible",
69  field_name_x.c_str(), field_name_y.c_str());
70  }
71 
72  typedef multi_index_container<
73  boost::shared_ptr<DofEntity>,
74  indexed_by<
75 
76  hashed_non_unique<
77  tag<Composite_Ent_Order_And_CoeffIdx_mi_tag>,
78  composite_key<
79 
80  DofEntity,
81  const_mem_fun<DofEntity, EntityHandle, &DofEntity::getEnt>,
82  const_mem_fun<DofEntity, ApproximationOrder,
84  const_mem_fun<DofEntity, FieldCoefficientsNumber,
86 
87  >>
88 
89  >>
90  DofEntity_multiIndex_composite_view;
91 
92  auto dof_lo_for_view =
93  dofs_ptr->get<FieldName_mi_tag>().lower_bound(field_name_y);
94  auto dof_hi_for_view =
95  dofs_ptr->get<FieldName_mi_tag>().upper_bound(field_name_y);
96  DofEntity_multiIndex_composite_view dof_composite_view;
97  dof_composite_view.insert(dof_lo_for_view, dof_hi_for_view);
98 
99  auto x_eit = field_ents->get<FieldName_mi_tag>().lower_bound(field_name_x);
100  auto x_eit_hi = field_ents->get<FieldName_mi_tag>().upper_bound(field_name_x);
101  for (; x_eit != x_eit_hi; x_eit++) {
102  int nb_dofs_on_x_entity = (*x_eit)->getNbDofsOnEnt();
103  VectorAdaptor field_data = (*x_eit)->getEntFieldData();
104  for (int dd = 0; dd != nb_dofs_on_x_entity; ++dd) {
105  ApproximationOrder dof_order = (*x_eit)->getDofOrderMap()[dd];
106  FieldCoefficientsNumber dof_rank = dd % (*x_eit)->getNbOfCoeffs();
107  FieldData data = field_data[dd];
108  auto dit = dof_composite_view.find(
109  boost::make_tuple((*x_eit)->getEnt(), dof_order, dof_rank));
110  if (dit == dof_composite_view.end()) {
111  if (creat_if_missing) {
112  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
113  "not yet implemented");
114  } else {
115  if (error_if_missing) {
116  std::ostringstream ss;
117  ss << "dof on ent " << (*x_eit)->getEnt() << " order " << dof_order
118  << " rank " << dof_rank << " does not exist";
119  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
120  ss.str().c_str());
121  } else {
122  continue;
123  }
124  }
125  }
126  CHKERR lambda((*dit)->getFieldData(),data);
127  }
128  }
130 }
131 
133  const std::string field_name_x,
134  const std::string field_name_y,
135  bool error_if_missing,
136  bool creat_if_missing) {
138  struct Axpy {
139  const double alpha;
140  Axpy(const double alpha) : alpha(alpha) {}
141  inline MoFEMErrorCode operator()(double &fy, const double fx) {
143  fy += alpha * fx;
145  }
146  };
147  CHKERR fieldLambda(Axpy(alpha), field_name_x, field_name_y, error_if_missing,
148  creat_if_missing);
150 }
151 
153  const std::string field_name_x,
154  const std::string field_name_y,
155  bool error_if_missing,
156  bool creat_if_missing) {
158  struct Copy {
159  const double alpha;
160  Copy(const double alpha) : alpha(alpha) {}
161  inline MoFEMErrorCode operator()(double &fy, const double fx) {
163  fy = alpha * fx;
165  }
166  };
167  CHKERR fieldLambda(Copy(alpha), field_name_x, field_name_y, error_if_missing,
168  creat_if_missing);
170 }
171 
173  const std::string field_name,
174  Range *sub_verts) {
175  const MoFEM::Interface &m_field = cOre;
177 
178  EntityHandle meshset = m_field.get_field_meshset(field_name);
179  Range verts;
180  CHKERR m_field.get_moab().get_entities_by_type(meshset, MBVERTEX, verts,
181  true);
182  if (sub_verts)
183  verts = intersect(*sub_verts, verts);
184 
185  struct LambdaMethod : EntityMethod {
186  LambdaMethod(MoFEM::Interface &m_field, Range &verts,
188  : EntityMethod(), mField(m_field), verts(verts), lambda(lambda),
189  count(0), total(0) {}
190  MoFEMErrorCode preProcess() {
191  vit = verts.begin();
192  return 0;
193  }
194  MoFEMErrorCode operator()() {
196  if (*vit != entPtr->getEnt())
197  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
198  "Inconsistent entity %ld != %ld", *vit, entPtr->getEnt());
199  if(!count)
200  CHKERR mField.get_moab().coords_iterate(vit, verts.end(), xCoord,
201  yCoord, zCoord, count);
202  CHKERR lambda(entPtr->getEntFieldData(), xCoord, yCoord, zCoord);
203  ++xCoord;
204  ++yCoord;
205  ++zCoord;
206  ++vit;
207  ++total;
208  --count;
210  }
211  MoFEMErrorCode postProcess() {
213  if(total != verts.size())
214  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
215  "Inconsistent number of vertices in field meshset and in the "
216  "field %d != %d",
217  total, verts.size());
219  }
220 
221  private:
222  MoFEM::Interface &mField;
223  Range &verts;
225  int count;
226  int total;
227  Range::iterator vit;
228  double *xCoord;
229  double *yCoord;
230  double *zCoord;
231  };
232 
233  LambdaMethod lambda_method(const_cast<MoFEM::Interface &>(m_field), verts,
234  lambda);
235  CHKERR const_cast<MoFEM::Interface &>(m_field).loop_entities(
236  field_name, lambda_method, &verts, QUIET);
238 }
239 
240 MoFEMErrorCode FieldBlas::setField(const double val, const EntityType type,
241  const std::string field_name) {
242  const MoFEM::Interface &m_field = cOre;
243  auto dofs_ptr = m_field.get_dofs();
245 
246  DofEntity_multiIndex::index<Composite_Name_And_Type_mi_tag>::type::iterator
247  dit,
248  hi_dit;
249  dit = dofs_ptr->get<Composite_Name_And_Type_mi_tag>().lower_bound(
250  boost::make_tuple(field_name, type));
251  hi_dit = dofs_ptr->get<Composite_Name_And_Type_mi_tag>().upper_bound(
252  boost::make_tuple(field_name, type));
253  for (; dit != hi_dit; dit++) {
254  (*dit)->getFieldData() = val;
255  }
257 }
258 
259 MoFEMErrorCode FieldBlas::setField(const double val, const EntityType type,
260  const Range &ents,
261  const std::string field_name) {
262  const MoFEM::Interface &m_field = cOre;
263  auto dofs_ptr = m_field.get_dofs();
265 
266  DofEntity_multiIndex::index<Composite_Name_And_Type_mi_tag>::type::iterator
267  dit,
268  hi_dit;
269  dit = dofs_ptr->get<Composite_Name_And_Type_mi_tag>().lower_bound(
270  boost::make_tuple(field_name, type));
271  hi_dit = dofs_ptr->get<Composite_Name_And_Type_mi_tag>().upper_bound(
272  boost::make_tuple(field_name, type));
273  EntityHandle ent, last = 0;
274  bool cont = true;
275  for (; dit != hi_dit; dit++) {
276  ent = (*dit)->getEnt();
277  if (ent != last) {
278  if (ents.find(ent) == ents.end()) {
279  cont = true;
280  } else {
281  cont = false;
282  }
283  last = ent;
284  }
285  if (cont)
286  continue;
287  (*dit)->getFieldData() = val;
288  }
290 }
291 
293  const std::string field_name) {
294  const MoFEM::Interface &m_field = cOre;
295  auto fields_ptr = m_field.get_fields();
296  auto dofs_ptr = m_field.get_dofs();
298  auto fit = fields_ptr->get<FieldName_mi_tag>().find(field_name);
299  if (fit == fields_ptr->get<FieldName_mi_tag>().end()) {
300  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
301  " field < %s > not found, (top tip: check spelling)",
302  field_name.c_str());
303  }
304 
305  auto dit = dofs_ptr->get<FieldName_mi_tag>().lower_bound(field_name);
306  auto hi_dit = dofs_ptr->get<FieldName_mi_tag>().upper_bound(field_name);
307  for (; dit != hi_dit; dit++) {
308  (*dit)->getFieldData() = val;
309  }
311 }
312 
314  const std::string field_name) {
315  const MoFEM::Interface &m_field = cOre;
316  auto fields_ptr = m_field.get_fields();
317  auto dofs_ptr = m_field.get_dofs();
319 
320  auto fit = fields_ptr->get<FieldName_mi_tag>().find(field_name);
321  if (fit == fields_ptr->get<FieldName_mi_tag>().end()) {
322  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
323  " field < %s > not found, (top tip: check spelling)",
324  field_name.c_str());
325  }
326 
327  DofEntityByFieldName::iterator dit, hi_dit;
328  dit = dofs_ptr->get<FieldName_mi_tag>().lower_bound(field_name);
329  hi_dit = dofs_ptr->get<FieldName_mi_tag>().upper_bound(field_name);
330  for (; dit != hi_dit; dit++) {
331  (*dit)->getFieldData() *= alpha;
332  }
334 }
335 
336 } // namespace MoFEM
Deprecated interface functions.
MoFEM interface unique ID.
FieldBlas(const MoFEM::Core &core)
Definition: FieldBlas.cpp:34
MoFEMErrorCode setField(const double val, const EntityType type, const std::string field_name)
scale field
Definition: FieldBlas.cpp:240
virtual moab::Interface & get_moab()=0
boost::function< MoFEMErrorCode(double &, const double)> TwoFieldFunction
Definition: FieldBlas.hpp:50
FieldCoefficientsNumber getDofCoeffIdx() const
virtual const DofEntity_multiIndex * get_dofs() const =0
Get the dofs object.
keeps information about DOF on the entity
int ApproximationOrder
Approximation on the entity.
Definition: Types.hpp:37
Data structure to exchange data between mofem and User Loop Methods on entities.It allows to exchange...
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:507
virtual const Field_multiIndex * get_fields() const =0
Get the fields object.
virtual const FieldEntity_multiIndex * get_field_ents() const =0
Get the field ents object.
base class for all interface classes
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
Core (interface) class.
Definition: Core.hpp:50
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:514
boost::function< MoFEMErrorCode(VectorAdaptor &&field_data, double *xcoord, double *ycoord, double *zcoord)> VertexCoordsFunction
Definition: FieldBlas.hpp:124
MoFEMErrorCode fieldAxpy(const double alpha, const std::string field_name_x, const std::string field_name_y, bool error_if_missing=false, bool creat_if_missing=false)
axpy fieldsfield_y = field_y + alpha*field_x
Definition: FieldBlas.cpp:132
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
MoFEMErrorCode query_interface(const MOFEMuuid &uuid, UnknownInterface **iface) const
Definition: FieldBlas.cpp:22
const MoFEM::Interface & cOre
Definition: FieldBlas.hpp:39
MoFEMErrorCode fieldScale(const double alpha, const std::string field_name)
set fieldfield_y = val
Definition: FieldBlas.cpp:313
~FieldBlas()
Destructor.
Definition: FieldBlas.cpp:36
virtual EntityHandle get_field_meshset(const std::string &name) const =0
get field meshset
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
MoFEMErrorCode fieldLambda(TwoFieldFunction lambda, const std::string field_name_x, const std::string field_name_y, bool error_if_missing=false, bool creat_if_missing=false)
filed lambdaDo calculation on two fields and save result to field fy
Definition: FieldBlas.cpp:38
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:108
#define CHKERR
Inline error check.
Definition: definitions.h:602
MultiIndex Tag for field name.
ApproximationOrder getDofOrder() const
static const MOFEMuuid IDD_MOFEMFieldBlas
Definition: FieldBlas.hpp:26
int FieldCoefficientsNumber
Number of field coefficients.
Definition: Types.hpp:38
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413
MoFEMErrorCode fieldCopy(const double alpha, const std::string field_name_x, const std::string field_name_y, bool error_if_missing=false, bool creat_if_missing=false)
copy and scale fieldsfield_y = alpha*field_x
Definition: FieldBlas.cpp:152
MoFEMErrorCode setVertexDofs(VertexCoordsFunction lambda, const std::string field_name, Range *verts=nullptr)
Set DOFs on vertices using user functionExample:
Definition: FieldBlas.cpp:172