v0.9.1
FieldApproximation.hpp
Go to the documentation of this file.
1 /* \file FieldApproximation.hpp
2 
3 \brief Element to calculate approximation on volume elements
4 
5 */
6 
7 /* This file is part of MoFEM.
8  * MoFEM is free software: you can redistribute it and/or modify it under
9  * the terms of the GNU Lesser General Public License as published by the
10  * Free Software Foundation, either version 3 of the License, or (at your
11  * option) any later version.
12  *
13  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
14  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  * License for more details.
17  *
18  * You should have received a copy of the GNU Lesser General Public
19  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
20 
21 #ifndef __FILEDAPPROXIMATION_HPP
22 #define __FILEDAPPROXIMATION_HPP
23 
24 using namespace boost::numeric;
25 
26 /** \brief Finite element for approximating analytical filed on the mesh
27  * \ingroup user_modules
28  */
30 
32  const std::string problemName;
35 
37  : mField(m_field), feVolume(m_field), feFace(m_field) {}
38 
39  /** \brief Gauss point operators to calculate matrices and vectors
40  *
41  * Function work on volumes (Terahedrons & Bricks)
42  */
43  template <typename FUNEVAL>
46 
47  Mat A;
48  std::vector<Vec> &vecF;
50 
51  OpApproxVolume(const std::string &field_name, Mat _A,
52  std::vector<Vec> &vec_F, FUNEVAL &function_evaluator)
54  field_name, UserDataOperator::OPROW | UserDataOperator::OPROWCOL),
55  A(_A), vecF(vec_F), functionEvaluator(function_evaluator) {}
56  virtual ~OpApproxVolume() {}
57 
60  std::vector<VectorDouble> Nf;
61 
62  /** \brief calculate matrix
63  */
64  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
65  EntityType col_type,
69 
70  if (A == PETSC_NULL)
72  if (row_data.getIndices().size() == 0)
74  if (col_data.getIndices().size() == 0)
76 
77  const auto &dof_ptr = row_data.getFieldDofs()[0];
78  int rank = dof_ptr->getNbOfCoeffs();
79 
80  int nb_row_dofs = row_data.getIndices().size() / rank;
81  int nb_col_dofs = col_data.getIndices().size() / rank;
82 
83  NN.resize(nb_row_dofs, nb_col_dofs, false);
84  NN.clear();
85 
86  unsigned int nb_gauss_pts = row_data.getN().size1();
87  for (unsigned int gg = 0; gg != nb_gauss_pts; gg++) {
88  double w = getVolume() * getGaussPts()(3, gg);
89  if (getHoCoordsAtGaussPts().size1() == nb_gauss_pts) {
90  w *= getHoGaussPtsDetJac()[gg];
91  }
92  // noalias(NN) += w*outer_prod(row_data.getN(gg),col_data.getN(gg));
93  cblas_dger(CblasRowMajor, nb_row_dofs, nb_col_dofs, w,
94  &row_data.getN()(gg, 0), 1, &col_data.getN()(gg, 0), 1,
95  &*NN.data().begin(), nb_col_dofs);
96  }
97 
98  if ((row_type != col_type) || (row_side != col_side)) {
99  transNN.resize(nb_col_dofs, nb_row_dofs, false);
100  ublas::noalias(transNN) = trans(NN);
101  }
102 
103  double *data = &*NN.data().begin();
104  double *trans_data = &*transNN.data().begin();
105  VectorInt row_indices, col_indices;
106  row_indices.resize(nb_row_dofs);
107  col_indices.resize(nb_col_dofs);
108 
109  for (int rr = 0; rr < rank; rr++) {
110 
111  if ((row_data.getIndices().size() % rank) != 0) {
112  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
113  "data inconsistency");
114  }
115 
116  if ((col_data.getIndices().size() % rank) != 0) {
117  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
118  "data inconsistency");
119  }
120 
121  unsigned int nb_rows;
122  unsigned int nb_cols;
123  int *rows;
124  int *cols;
125 
126  if (rank > 1) {
127 
128  ublas::noalias(row_indices) = ublas::vector_slice<VectorInt>(
129  row_data.getIndices(),
130  ublas::slice(rr, rank, row_data.getIndices().size() / rank));
131  ublas::noalias(col_indices) = ublas::vector_slice<VectorInt>(
132  col_data.getIndices(),
133  ublas::slice(rr, rank, col_data.getIndices().size() / rank));
134 
135  nb_rows = row_indices.size();
136  nb_cols = col_indices.size();
137  rows = &*row_indices.data().begin();
138  cols = &*col_indices.data().begin();
139 
140  } else {
141 
142  nb_rows = row_data.getIndices().size();
143  nb_cols = col_data.getIndices().size();
144  rows = &*row_data.getIndices().data().begin();
145  cols = &*col_data.getIndices().data().begin();
146  }
147 
148  if (nb_rows != NN.size1()) {
149  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
150  "data inconsistency");
151  }
152  if (nb_cols != NN.size2()) {
153  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
154  "data inconsistency");
155  }
156 
157  CHKERR MatSetValues(A, nb_rows, rows, nb_cols, cols, data, ADD_VALUES);
158  if ((row_type != col_type) || (row_side != col_side)) {
159  if (nb_rows != transNN.size2()) {
160  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
161  "data inconsistency");
162  }
163  if (nb_cols != transNN.size1()) {
164  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
165  "data inconsistency");
166  }
167  CHKERR MatSetValues(A, nb_cols, cols, nb_rows, rows, trans_data,
168  ADD_VALUES);
169  }
170  }
171 
173  }
174 
175  /** \brief calculate vector
176  */
177  MoFEMErrorCode doWork(int side, EntityType type,
180 
181  if (data.getIndices().size() == 0)
183 
184  // PetscAttachDebugger();
185 
186  const auto &dof_ptr = data.getFieldDofs()[0];
187  unsigned int rank = dof_ptr->getNbOfCoeffs();
188 
189  int nb_row_dofs = data.getIndices().size() / rank;
190 
191  if (getCoordsAtGaussPts().size1() != data.getN().size1()) {
192  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
193  "data inconsistency");
194  }
195  if (getCoordsAtGaussPts().size2() != 3) {
196  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
197  "data inconsistency");
198  }
199 
200  // itegration
201  unsigned int nb_gauss_pts = data.getN().size1();
202  for (unsigned int gg = 0; gg != nb_gauss_pts; gg++) {
203 
204  double x, y, z, w;
205  w = getVolume() * getGaussPts()(3, gg);
206  if (getHoCoordsAtGaussPts().size1() == nb_gauss_pts) {
207  // intergation points global positions if higher order geometry is
208  // given
209  x = getHoCoordsAtGaussPts()(gg, 0);
210  y = getHoCoordsAtGaussPts()(gg, 1);
211  z = getHoCoordsAtGaussPts()(gg, 2);
212  // correction of jacobian for higher order geometry
213  w *= getHoGaussPtsDetJac()[gg];
214  } else {
215  // intergartion point global positions for linear tetrahedral element
216  x = getCoordsAtGaussPts()(gg, 0);
217  y = getCoordsAtGaussPts()(gg, 1);
218  z = getCoordsAtGaussPts()(gg, 2);
219  }
220 
221  std::vector<VectorDouble> fun_val;
222 
223  fun_val = functionEvaluator(x, y, z);
224 
225  if (fun_val.size() != vecF.size()) {
226  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
227  "data inconsistency");
228  }
229 
230  Nf.resize(fun_val.size());
231  for (unsigned int lhs = 0; lhs != fun_val.size(); lhs++) {
232 
233  if (!gg) {
234  Nf[lhs].resize(data.getIndices().size());
235  Nf[lhs].clear();
236  }
237 
238  if (fun_val[lhs].size() != rank) {
239  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
240  "data inconsistency");
241  }
242 
243  for (unsigned int rr = 0; rr != rank; rr++) {
244  cblas_daxpy(nb_row_dofs, w * (fun_val[lhs])[rr],
245  &data.getN()(gg, 0), 1, &(Nf[lhs])[rr], rank);
246  }
247  }
248  }
249 
250  for (unsigned int lhs = 0; lhs != vecF.size(); lhs++) {
251 
252  CHKERR VecSetValues(vecF[lhs], data.getIndices().size(),
253  &data.getIndices()[0], &(Nf[lhs])[0], ADD_VALUES);
254  }
255 
257  }
258  };
259 
260  /** \brief Gauss point operators to calculate matrices and vectors
261  *
262  * Function work on faces (Triangles & Quads)
263  */
264  template <typename FUNEVAL>
267 
268  Mat A;
269  std::vector<Vec> &vecF;
271 
272  OpApproxFace(const std::string &field_name, Mat _A, std::vector<Vec> &vec_F,
273  FUNEVAL &function_evaluator)
275  field_name, UserDataOperator::OPROW | UserDataOperator::OPROWCOL),
276  A(_A), vecF(vec_F), functionEvaluator(function_evaluator) {}
277  virtual ~OpApproxFace() {}
278 
281  std::vector<VectorDouble> Nf;
282 
283  /** \brief calculate matrix
284  */
285  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
286  EntityType col_type,
290 
291  if (A == PETSC_NULL)
293  if (row_data.getIndices().size() == 0)
295  if (col_data.getIndices().size() == 0)
297 
298  const auto &dof_ptr = row_data.getFieldDofs()[0];
299  int rank = dof_ptr->getNbOfCoeffs();
300  int nb_row_dofs = row_data.getIndices().size() / rank;
301  int nb_col_dofs = col_data.getIndices().size() / rank;
302  NN.resize(nb_row_dofs, nb_col_dofs, false);
303  NN.clear();
304  unsigned int nb_gauss_pts = row_data.getN().size1();
305  for (unsigned int gg = 0; gg != nb_gauss_pts; gg++) {
306  double w = getGaussPts()(2, gg);
307  if (getNormalsAtGaussPts().size1()) {
308  w *= 0.5 * cblas_dnrm2(3, &getNormalsAtGaussPts()(gg, 0), 1);
309  } else {
310  w *= getArea();
311  }
312  // noalias(NN) += w*outer_prod(row_data.getN(gg),col_data.getN(gg));
313  cblas_dger(CblasRowMajor, nb_row_dofs, nb_col_dofs, w,
314  &row_data.getN()(gg, 0), 1, &col_data.getN()(gg, 0), 1,
315  &*NN.data().begin(), nb_col_dofs);
316  }
317 
318  if ((row_type != col_type) || (row_side != col_side)) {
319  transNN.resize(nb_col_dofs, nb_row_dofs, false);
320  ublas::noalias(transNN) = trans(NN);
321  }
322 
323  double *data = &*NN.data().begin();
324  double *trans_data = &*transNN.data().begin();
325  VectorInt row_indices, col_indices;
326  row_indices.resize(nb_row_dofs);
327  col_indices.resize(nb_col_dofs);
328  for (int rr = 0; rr < rank; rr++) {
329  if ((row_data.getIndices().size() % rank) != 0) {
330  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
331  "data inconsistency");
332  }
333  if ((col_data.getIndices().size() % rank) != 0) {
334  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
335  "data inconsistency");
336  }
337  unsigned int nb_rows;
338  unsigned int nb_cols;
339  int *rows;
340  int *cols;
341  if (rank > 1) {
342  ublas::noalias(row_indices) = ublas::vector_slice<VectorInt>(
343  row_data.getIndices(),
344  ublas::slice(rr, rank, row_data.getIndices().size() / rank));
345  ublas::noalias(col_indices) = ublas::vector_slice<VectorInt>(
346  col_data.getIndices(),
347  ublas::slice(rr, rank, col_data.getIndices().size() / rank));
348  nb_rows = row_indices.size();
349  nb_cols = col_indices.size();
350  rows = &*row_indices.data().begin();
351  cols = &*col_indices.data().begin();
352  } else {
353  nb_rows = row_data.getIndices().size();
354  nb_cols = col_data.getIndices().size();
355  rows = &*row_data.getIndices().data().begin();
356  cols = &*col_data.getIndices().data().begin();
357  }
358  if (nb_rows != NN.size1()) {
359  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
360  "data inconsistency");
361  }
362  if (nb_cols != NN.size2()) {
363  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
364  "data inconsistency");
365  }
366  CHKERR MatSetValues(A, nb_rows, rows, nb_cols, cols, data, ADD_VALUES);
367  if ((row_type != col_type) || (row_side != col_side)) {
368  if (nb_rows != transNN.size2()) {
369  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
370  "data inconsistency");
371  }
372  if (nb_cols != transNN.size1()) {
373  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
374  "data inconsistency");
375  }
376  CHKERR MatSetValues(A, nb_cols, cols, nb_rows, rows, trans_data,
377  ADD_VALUES);
378  }
379  }
381  }
382 
383  /** \brief calculate vector
384  */
385  MoFEMErrorCode doWork(int side, EntityType type,
388 
389  if (data.getIndices().size() == 0)
391 
392  const auto &dof_ptr = data.getFieldDofs()[0];
393  unsigned int rank = dof_ptr->getNbOfCoeffs();
394 
395  int nb_row_dofs = data.getIndices().size() / rank;
396 
397  if (getCoordsAtGaussPts().size1() != data.getN().size1()) {
398  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
399  "data inconsistency");
400  }
401  if (getCoordsAtGaussPts().size2() != 3) {
402  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
403  "data inconsistency");
404  }
405 
406  VectorDouble normal(3);
407  VectorDouble tangent1(3);
408  VectorDouble tangent2(3);
409  tangent1.clear();
410  tangent2.clear();
411 
412  // integration
413  unsigned int nb_gauss_pts = data.getN().size1();
414  for (unsigned int gg = 0; gg != nb_gauss_pts; gg++) {
415  double x, y, z, w;
416  w = getGaussPts()(2, gg);
417  if (getNormalsAtGaussPts().size1()) {
418  w *= 0.5 * cblas_dnrm2(3, &getNormalsAtGaussPts()(gg, 0), 1);
419  } else {
420  w *= getArea();
421  }
422 
423  if (getHoCoordsAtGaussPts().size1() == nb_gauss_pts) {
424  // intergation points global positions if higher order geometry is
425  // given
426  x = getHoCoordsAtGaussPts()(gg, 0);
427  y = getHoCoordsAtGaussPts()(gg, 1);
428  z = getHoCoordsAtGaussPts()(gg, 2);
429  } else {
430  // intergartion point global positions for linear tetrahedral element
431  x = getCoordsAtGaussPts()(gg, 0);
432  y = getCoordsAtGaussPts()(gg, 1);
433  z = getCoordsAtGaussPts()(gg, 2);
434  }
435 
436  if (getNormalsAtGaussPts().size1()) {
437  noalias(normal) = getNormalsAtGaussPts(gg);
438  for (int dd = 0; dd < 3; dd++) {
439  tangent1[dd] = getTangent1AtGaussPts()(gg, dd);
440  tangent2[dd] = getTangent2AtGaussPts()(gg, dd);
441  }
442  } else {
443  noalias(normal) = getNormal();
444  }
445 
446  std::vector<VectorDouble> fun_val;
447  EntityHandle ent = getFEMethod()->numeredEntFiniteElementPtr->getEnt();
448  fun_val = functionEvaluator(ent, x, y, z, normal, tangent1, tangent2);
449  if (fun_val.size() != vecF.size()) {
450  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
451  "data inconsistency");
452  }
453 
454  Nf.resize(fun_val.size());
455  for (unsigned int lhs = 0; lhs != fun_val.size(); lhs++) {
456  if (!gg) {
457  Nf[lhs].resize(data.getIndices().size());
458  Nf[lhs].clear();
459  }
460  if (fun_val[lhs].size() != rank) {
461  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
462  "data inconsistency");
463  }
464  for (unsigned int rr = 0; rr != rank; rr++) {
465  cblas_daxpy(nb_row_dofs, w * (fun_val[lhs])[rr],
466  &data.getN()(gg, 0), 1, &(Nf[lhs])[rr], rank);
467  }
468  }
469  }
470 
471  for (unsigned int lhs = 0; lhs != vecF.size(); lhs++) {
472 
473  CHKERR VecSetValues(vecF[lhs], data.getIndices().size(),
474  &data.getIndices()[0], &(Nf[lhs])[0], ADD_VALUES);
475  }
476 
478  }
479  };
480 
481  /** \brief Set operators
482  */
483  template <typename FUNEVAL>
484  MoFEMErrorCode setOperatorsVolume(const std::string &field_name, Mat A,
485  std::vector<Vec> &vec_F,
486  FUNEVAL &function_evaluator) {
488  // add operator to calculate F vector
489  feVolume.getOpPtrVector().push_back(
490  new OpApproxVolume<FUNEVAL>(field_name, A, vec_F, function_evaluator));
491  // add operator to calculate A matrix
492  // if(A) {
493  // feVolume.getOpPtrVector().push_back(new
494  // OpApproxVolume<FUNEVAL>(field_name,A,vec_F,function_evaluator));
495  // }
497  }
498 
499  /** \brief Set operators
500  */
501  template <typename FUNEVAL>
502  MoFEMErrorCode setOperatorsFace(const std::string &field_name, Mat A,
503  std::vector<Vec> &vec_F,
504  FUNEVAL &function_evaluator) {
506  // add operator to calculate F vector
507  feFace.getOpPtrVector().push_back(
508  new OpApproxFace<FUNEVAL>(field_name, A, vec_F, function_evaluator));
509  // add operator to calculate A matrix
510  // if(A) {
511  // feFace.getOpPtrVector().push_back(new
512  // OpApproxFace<FUNEVAL>(field_name,A,vec_F,function_evaluator));
513  // }
515  }
516 
517  /** \brief assemble matrix and vector
518  */
519  template <typename FUNEVAL>
520  MoFEMErrorCode loopMatrixAndVectorVolume(const std::string &problem_name,
521  const std::string &fe_name,
522  const std::string &field_name, Mat A,
523  std::vector<Vec> &vec_F,
524  FUNEVAL &function_evaluator) {
526 
527  CHKERR setOperatorsVolume(field_name, A, vec_F, function_evaluator);
528  if (A) {
529  CHKERR MatZeroEntries(A);
530  }
531  // calculate and assemble
532  CHKERR mField.loop_finite_elements(problem_name, fe_name, feVolume);
533  if (A) {
534  CHKERR MatAssemblyBegin(A, MAT_FLUSH_ASSEMBLY);
535  CHKERR MatAssemblyEnd(A, MAT_FLUSH_ASSEMBLY);
536  }
537  for (unsigned int lhs = 0; lhs < vec_F.size(); lhs++) {
538  CHKERR VecAssemblyBegin(vec_F[lhs]);
539  CHKERR VecAssemblyEnd(vec_F[lhs]);
540  }
542  }
543 };
544 
545 #endif //__FILEDAPPROXIMATION_HPP
Deprecated interface functions.
void cblas_daxpy(const int N, const double alpha, const double *X, const int incX, double *Y, const int incY)
Definition: cblas_daxpy.c:11
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
calculate matrix
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
FTensor::Tensor1< double, 2 > normal(FTensor::Tensor1< T1, 3 > &t_coords, FTensor::Tensor1< T2, 2 > &t_disp)
Definition: ContactOps.hpp:128
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
void cblas_dger(const enum CBLAS_ORDER order, const int M, const int N, const double alpha, const double *X, const int incX, const double *Y, const int incY, double *A, const int lda)
Definition: cblas_dger.c:12
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:513
MoFEMErrorCode setOperatorsVolume(const std::string &field_name, Mat A, std::vector< Vec > &vec_F, FUNEVAL &function_evaluator)
Set operators.
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
calculate vector
double cblas_dnrm2(const int N, const double *X, const int incX)
Definition: cblas_dnrm2.c:12
MoFEM::Interface & mField
Gauss point operators to calculate matrices and vectors.
MoFEMErrorCode loopMatrixAndVectorVolume(const std::string &problem_name, const std::string &fe_name, const std::string &field_name, Mat A, std::vector< Vec > &vec_F, FUNEVAL &function_evaluator)
assemble matrix and vector
FieldApproximationH1(MoFEM::Interface &m_field)
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
MoFEM::FaceElementForcesAndSourcesCore feFace
MoFEMErrorCode setOperatorsFace(const std::string &field_name, Mat A, std::vector< Vec > &vec_F, FUNEVAL &function_evaluator)
Set operators.
OpApproxFace(const std::string &field_name, Mat _A, std::vector< Vec > &vec_F, FUNEVAL &function_evaluator)
DataForcesAndSourcesCore::EntData EntData
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
MoFEMErrorCode MatSetValues(Mat M, const DataForcesAndSourcesCore::EntData &row_data, const DataForcesAndSourcesCore::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
MoFEMErrorCode VecSetValues(Vec V, const DataForcesAndSourcesCore::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
calculate matrix
#define CHKERR
Inline error check.
Definition: definitions.h:601
ForcesAndSourcesCore::UserDataOperator UserDataOperator
VolumeElementForcesAndSourcesCore feVolume
ublas::vector< int, IntAllocator > VectorInt
Definition: Types.hpp:71
const std::string problemName
Finite element for approximating analytical filed on the mesh.
Gauss point operators to calculate matrices and vectors.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
calculate vector
OpApproxVolume(const std::string &field_name, Mat _A, std::vector< Vec > &vec_F, FUNEVAL &function_evaluator)
DEPRECATED MoFEMErrorCode loop_finite_elements(const Problem *problem_ptr, const std::string &fe_name, FEMethod &method, int lower_rank, int upper_rank, MoFEMTypes bh, int verb=DEFAULT_VERBOSITY)
double w(const double g, const double t)
Definition: ContactOps.hpp:160