v0.13.0
scalar_check_approximation.cpp
Go to the documentation of this file.
1 /**
2  * \file scalar_check_approximation.cpp
3  * \example scalar_check_approximation.cpp
4  *
5  * Approximate polynomial in 2D and check derivatives
6  *
7  */
8 
9 /* This file is part of MoFEM.
10  * MoFEM is free software: you can redistribute it and/or modify it under
11  * the terms of the GNU Lesser General Public License as published by the
12  * Free Software Foundation, either version 3 of the License, or (at your
13  * option) any later version.
14  *
15  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
16  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Publicrule
18  * License for more details.
19  *
20  * You should have received a copy of the GNU Lesser General Public
21  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
22 
23 #include <MoFEM.hpp>
24 
25 using namespace MoFEM;
26 
27 static char help[] = "...\n\n";
28 
29 static int approx_order = 4;
30 
31 template <int DIM> struct ApproxFunctionsImpl {};
32 
33 template <int DIM> struct ElementsAndOps {};
34 
35 template <> struct ElementsAndOps<2> {
38 };
39 
40 template <> struct ElementsAndOps<3> {
43 };
44 
45 constexpr int SPACE_DIM =
46  EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
47 
51 
52 template <> struct ApproxFunctionsImpl<2> {
53  static double fUn(const double x, const double y, double z) {
54  double r = 1;
55  for (int o = 1; o <= approx_order; ++o) {
56  for (int i = 0; i <= o; ++i) {
57  int j = o - i;
58  if (j >= 0)
59  r += pow(x, i) * pow(y, j);
60  }
61  }
62  return r;
63  }
64 
65  static FTensor::Tensor1<double, 2> diffFun(const double x, const double y,
66  double z) {
68  for (int o = 1; o <= approx_order; ++o) {
69  for (int i = 0; i <= o; ++i) {
70  int j = o - i;
71  if (j >= 0) {
72  r(0) += i > 0 ? i * pow(x, i - 1) * pow(y, j) : 0;
73  r(1) += j > 0 ? j * pow(x, i) * pow(y, j - 1) : 0;
74  }
75  }
76  }
77  return r;
78  }
79 };
80 
81 template <> struct ApproxFunctionsImpl<3> {
82  static double fUn(const double x, const double y, double z) {
83  double r = 1;
84  for (int o = 1; o <= approx_order; ++o) {
85  for (int i = 0; i <= o; ++i) {
86  for (int j = 0; j <= o - i; j++) {
87  int k = o - i - j;
88  if (k >= 0) {
89  r += pow(x, i) * pow(y, j) * pow(z, k);
90  }
91  }
92  }
93  }
94  return r;
95  }
96 
97  static FTensor::Tensor1<double, 3> diffFun(const double x, const double y,
98  double z) {
99  FTensor::Tensor1<double, 3> r{0., 0., 0.};
100  for (int o = 1; o <= approx_order; ++o) {
101  for (int i = 0; i <= o; ++i) {
102  for (int j = 0; j <= o - i; j++) {
103  int k = o - i - j;
104  if (k >= 0) {
105  r(0) += i > 0 ? i * pow(x, i - 1) * pow(y, j) * pow(z, k) : 0;
106  r(1) += j > 0 ? j * pow(x, i) * pow(y, j - 1) * pow(z, k) : 0;
107  r(2) += k > 0 ? k * pow(x, i) * pow(y, j) * pow(z, k - 1) : 0;
108  }
109  }
110  }
111  }
112  return r;
113  }
114 };
115 
117 
118 struct OpValsDiffVals : public DomainEleOp {
120  MatrixDouble &diffVals;
121  const bool checkGradients;
122  OpValsDiffVals(VectorDouble &vals, MatrixDouble &diff_vals, bool check_grads)
123  : DomainEleOp("FIELD1", OPROW), vAls(vals), diffVals(diff_vals),
124  checkGradients(check_grads) {}
125 
127 
128  MoFEMErrorCode doWork(int side, EntityType type,
131  const int nb_gauss_pts = getGaussPts().size2();
132  if (type == MBVERTEX) {
133  vAls.resize(nb_gauss_pts, false);
134  diffVals.resize(SPACE_DIM, nb_gauss_pts, false);
135  vAls.clear();
136  diffVals.clear();
137  }
138 
139  const int nb_dofs = data.getIndices().size();
140  if (nb_dofs) {
141 
142  MOFEM_LOG("AT", Sev::noisy)
143  << "Type " << moab::CN::EntityTypeName(type) << " side " << side;
144  MOFEM_LOG("AT", Sev::noisy) << data.getN();
145  MOFEM_LOG("AT", Sev::noisy) << data.getDiffN();
146 
147  auto t_vals = getFTensor0FromVec(vAls);
148  auto t_base_fun = data.getFTensor0N();
149  for (int gg = 0; gg != nb_gauss_pts; gg++) {
150  auto t_data = data.getFTensor0FieldData();
151  for (int bb = 0; bb != nb_dofs; bb++) {
152  t_vals += t_base_fun * t_data;
153  ++t_base_fun;
154  ++t_data;
155  }
156  const double v = t_vals;
157  if (!std::isnormal(v))
158  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "Not a number");
159  ++t_vals;
160  }
161 
162  if (checkGradients) {
163  auto t_diff_vals = getFTensor1FromMat<SPACE_DIM>(diffVals);
164  auto t_diff_base_fun = data.getFTensor1DiffN<SPACE_DIM>();
165  for (int gg = 0; gg != nb_gauss_pts; gg++) {
166  auto t_data = data.getFTensor0FieldData();
167  for (int bb = 0; bb != nb_dofs; bb++) {
168  t_diff_vals(i) += t_diff_base_fun(i) * t_data;
169  ++t_diff_base_fun;
170  ++t_data;
171  }
172  for (int d = 0; d != SPACE_DIM; ++d)
173  if (!std::isnormal(t_diff_vals(d)))
174  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "Not a number");
175  ++t_diff_vals;
176  }
177  }
178  }
180  }
181 };
182 
183 struct OpCheckValsDiffVals : public DomainEleOp {
185  MatrixDouble &diffVals;
186  boost::shared_ptr<VectorDouble> ptrVals;
187  boost::shared_ptr<MatrixDouble> ptrDiffVals;
188  const bool checkGradients;
189 
191  boost::shared_ptr<VectorDouble> &ptr_vals,
192  boost::shared_ptr<MatrixDouble> &ptr_diff_vals,
193  bool check_grads)
194  : DomainEleOp("FIELD1", OPROW), vAls(vals), diffVals(diff_vals),
195  ptrVals(ptr_vals), ptrDiffVals(ptr_diff_vals),
196  checkGradients(check_grads) {
197  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
198  }
199 
201 
202  MoFEMErrorCode doWork(int side, EntityType type,
205  const double eps = 1e-6;
206  const int nb_gauss_pts = data.getN().size1();
207 
208  auto t_vals = getFTensor0FromVec(vAls);
209 
210  auto t_ptr_vals = getFTensor0FromVec(*ptrVals);
211 
212  for (int gg = 0; gg != nb_gauss_pts; gg++) {
213 
214  double err_val;
215 
216  // Check user data operators
217  err_val = std::abs(t_vals - t_ptr_vals);
218  MOFEM_LOG("AT", Sev::noisy) << "Val op error " << err_val;
219 
220  if (err_val > eps)
221  SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
222  "Wrong value from operator %4.3e", err_val);
223 
224  const double x = getCoordsAtGaussPts()(gg, 0);
225  const double y = getCoordsAtGaussPts()(gg, 1);
226  const double z = getCoordsAtGaussPts()(gg, 2);
227 
228  // Check approximation
229  const double delta_val = t_vals - ApproxFunctions::fUn(x, y, z);
230 
231  err_val = std::fabs(delta_val * delta_val);
232  MOFEM_LOG("AT", Sev::verbose) << err_val << " : " << t_vals;
233  if (err_val > eps)
234  SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "Wrong value %4.3e",
235  err_val);
236 
237  ++t_vals;
238  ++t_ptr_vals;
239  }
240 
241  if (checkGradients) {
242 
243  auto t_diff_vals = getFTensor1FromMat<SPACE_DIM>(diffVals);
244  auto t_ptr_diff_vals = getFTensor1FromMat<SPACE_DIM>(*ptrDiffVals);
245 
246  for (int gg = 0; gg != nb_gauss_pts; gg++) {
247 
248  FTensor::Tensor1<double, SPACE_DIM> t_delta_diff_val;
249  double err_diff_val;
250 
251  t_delta_diff_val(i) = t_diff_vals(i) - t_ptr_diff_vals(i);
252  err_diff_val = sqrt(t_delta_diff_val(i) * t_delta_diff_val(i));
253  MOFEM_LOG("AT", Sev::noisy) << "Diff val op error " << err_diff_val;
254 
255  if (err_diff_val > eps)
256  SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
257  "Wrong derivatives from operator %4.3e", err_diff_val);
258 
259  const double x = getCoordsAtGaussPts()(gg, 0);
260  const double y = getCoordsAtGaussPts()(gg, 1);
261  const double z = getCoordsAtGaussPts()(gg, 2);
262 
263  // Check approximation
264  auto t_diff_anal = ApproxFunctions::diffFun(x, y, z);
265  t_delta_diff_val(i) = t_diff_vals(i) - t_diff_anal(i);
266 
267  err_diff_val = sqrt(t_delta_diff_val(i) * t_delta_diff_val(i));
268  if (SPACE_DIM == 3)
269  MOFEM_LOG("AT", Sev::noisy)
270  << "Diff val " << err_diff_val << " : "
271  << sqrt(t_diff_vals(i) * t_diff_vals(i)) << " : "
272  << t_diff_vals(0) << " (" << t_diff_anal(0) << ") "
273  << t_diff_vals(1) << " (" << t_diff_anal(1) << ") "
274  << t_diff_vals(2) << " (" << t_diff_anal(2) << ")";
275  else
276  MOFEM_LOG("AT", Sev::noisy)
277  << "Diff val " << err_diff_val << " : "
278  << sqrt(t_diff_vals(i) * t_diff_vals(i)) << " : "
279  << t_diff_vals(0) << " (" << t_diff_anal(0) << ") "
280  << t_diff_vals(1) << " (" << t_diff_anal(1) << ")";
281 
282  MOFEM_LOG("AT", Sev::verbose)
283  << getCoords()(3 * 1 + 0) - getCoords()(3 * 0 + 0);
284  MOFEM_LOG("AT", Sev::verbose)
285  << getCoords()(3 * 1 + 1) - getCoords()(3 * 0 + 1);
286  MOFEM_LOG("AT", Sev::verbose)
287  << getCoords()(3 * 1 + 2) - getCoords()(3 * 0 + 2);
288 
289  MOFEM_LOG("AT", Sev::verbose) << "Diff val error " << err_diff_val;
290  if (err_diff_val > eps)
291  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
292  "Wrong derivative of value %4.3e %4.3e", err_diff_val,
293  t_diff_anal.l2());
294 
295  ++t_diff_vals;
296  ++t_ptr_diff_vals;
297  }
298  }
299 
301  }
302 };
303 
304 int main(int argc, char *argv[]) {
305 
306  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
307 
308  try {
309 
310  DMType dm_name = "DMMOFEM";
311  CHKERR DMRegister_MoFEM(dm_name);
312 
313  moab::Core mb_instance;
314  moab::Interface &moab = mb_instance;
315 
316  // Add logging channel for example
317  auto core_log = logging::core::get();
318  core_log->add_sink(
320  LogManager::setLog("AT");
321  MOFEM_LOG_TAG("AT", "atom_test");
322 
323  // Create MoFEM instance
324  MoFEM::Core core(moab);
325  MoFEM::Interface &m_field = core;
326 
327  Simple *simple_interface = m_field.getInterface<Simple>();
328  PipelineManager *pipeline_mng = m_field.getInterface<PipelineManager>();
329  CHKERR simple_interface->getOptions();
330  CHKERR simple_interface->loadFile("", "");
331 
332  // Declare elements
333  enum bases {
334  AINSWORTH,
335  AINSWORTH_LOBATTO,
336  DEMKOWICZ,
337  BERNSTEIN,
338  LASBASETOP
339  };
340  const char *list_bases[] = {"ainsworth", "ainsworth_labatto", "demkowicz",
341  "bernstein"};
342  PetscBool flg;
343  PetscInt choice_base_value = AINSWORTH;
344  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
345  LASBASETOP, &choice_base_value, &flg);
346 
347  if (flg != PETSC_TRUE)
348  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "base not set");
350  if (choice_base_value == AINSWORTH)
352  if (choice_base_value == AINSWORTH_LOBATTO)
353  base = AINSWORTH_LOBATTO_BASE;
354  else if (choice_base_value == DEMKOWICZ)
355  base = DEMKOWICZ_JACOBI_BASE;
356  else if (choice_base_value == BERNSTEIN)
358 
359  enum spaces { H1SPACE, L2SPACE, LASBASETSPACE };
360  const char *list_spaces[] = {"h1", "l2"};
361  PetscInt choice_space_value = H1SPACE;
362  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-space", list_spaces,
363  LASBASETSPACE, &choice_space_value, &flg);
364  if (flg != PETSC_TRUE)
365  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "space not set");
366  FieldSpace space = H1;
367  if (choice_space_value == H1SPACE)
368  space = H1;
369  else if (choice_space_value == L2SPACE)
370  space = L2;
371 
372  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &approx_order,
373  PETSC_NULL);
374 
375  CHKERR simple_interface->addDomainField("FIELD1", space, base, 1);
376  CHKERR simple_interface->setFieldOrder("FIELD1", approx_order);
377 
378  Range edges, faces;
379  CHKERR moab.get_entities_by_dimension(0, 1, edges);
380  CHKERR moab.get_entities_by_dimension(0, 2, faces);
381 
382  if (choice_base_value != BERNSTEIN) {
383  Range rise_order;
384 
385  int i = 0;
386  for (auto e : edges) {
387  if (!(i % 2)) {
388  rise_order.insert(e);
389  }
390  ++i;
391  }
392 
393  for (auto f : faces) {
394  if (!(i % 3)) {
395  rise_order.insert(f);
396  }
397  ++i;
398  }
399 
400  Range rise_twice;
401  for (auto e : rise_order) {
402  if (!(i % 2)) {
403  rise_twice.insert(e);
404  }
405  ++i;
406  }
407 
408  CHKERR simple_interface->setFieldOrder("FIELD1", approx_order + 1,
409  &rise_order);
410 
411  CHKERR simple_interface->setFieldOrder("FIELD1", approx_order + 2,
412  &rise_twice);
413  }
414 
415  CHKERR simple_interface->setUp();
416  auto dm = simple_interface->getDM();
417 
418  VectorDouble vals;
419  auto jac_ptr = boost::make_shared<MatrixDouble>();
420  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
421  auto det_ptr = boost::make_shared<VectorDouble>();
422  MatrixDouble diff_vals;
423 
424  auto assemble_matrices_and_vectors = [&]() {
426  if (SPACE_DIM == 2)
427  pipeline_mng->getOpDomainRhsPipeline().push_back(
428  new OpSetHOWeightsOnFace());
429 
432 
433  if (SPACE_DIM == 2) {
434  pipeline_mng->getOpDomainLhsPipeline().push_back(
435  new OpCalculateHOJacForFace(jac_ptr));
436  pipeline_mng->getOpDomainLhsPipeline().push_back(
437  new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
438  pipeline_mng->getOpDomainLhsPipeline().push_back(
439  new OpSetHOWeightsOnFace());
440  }
441 
442  if (SPACE_DIM == 3) {
443  pipeline_mng->getOpDomainLhsPipeline().push_back(
444  new OpCalculateHOJacVolume(jac_ptr));
445  pipeline_mng->getOpDomainLhsPipeline().push_back(
446  new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, nullptr));
447  pipeline_mng->getOpDomainLhsPipeline().push_back(
448  new OpSetHOWeights(det_ptr));
449  pipeline_mng->getOpDomainRhsPipeline().push_back(
450  new OpCalculateHOJacVolume(jac_ptr));
451  pipeline_mng->getOpDomainRhsPipeline().push_back(
452  new OpInvertMatrix<3>(jac_ptr, det_ptr, nullptr));
453  pipeline_mng->getOpDomainRhsPipeline().push_back(
454  new OpSetHOWeights(det_ptr));
455  }
456 
459  pipeline_mng->getOpDomainLhsPipeline().push_back(new OpMass(
460  "FIELD1", "FIELD1", [](double, double, double) { return 1.; }));
461  pipeline_mng->getOpDomainRhsPipeline().push_back(
462  new OpSource("FIELD1", ApproxFunctions::fUn));
463 
464  auto integration_rule = [](int, int, int p_data) {
465  return 2 * p_data + 1;
466  };
467  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
468  CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
469 
471  };
472 
473  auto solve_problem = [&] {
475  auto solver = pipeline_mng->createKSP();
476  CHKERR KSPSetFromOptions(solver);
477  CHKERR KSPSetUp(solver);
478 
479  auto dm = simple_interface->getDM();
480  auto D = smartCreateDMVector(dm);
481  auto F = smartVectorDuplicate(D);
482 
483  CHKERR KSPSolve(solver, F, D);
484  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
485  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
486  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
487 
489  };
490 
491  auto check_solution = [&] {
493 
494  boost::shared_ptr<VectorDouble> ptr_values(new VectorDouble());
495  boost::shared_ptr<MatrixDouble> ptr_diff_vals(new MatrixDouble());
496 
497  pipeline_mng->getDomainLhsFE().reset();
498  pipeline_mng->getOpDomainRhsPipeline().clear();
499 
500  if (SPACE_DIM == 2) {
501  pipeline_mng->getOpDomainRhsPipeline().push_back(
502  new OpCalculateHOJacForFace(jac_ptr));
503  pipeline_mng->getOpDomainRhsPipeline().push_back(
504  new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
505  pipeline_mng->getOpDomainRhsPipeline().push_back(
506  new OpSetInvJacSpaceForFaceImpl<2>(space, inv_jac_ptr));
507  }
508 
509  if (SPACE_DIM == 3) {
510  pipeline_mng->getOpDomainRhsPipeline().push_back(
511  new OpCalculateHOJacVolume(jac_ptr));
512  pipeline_mng->getOpDomainRhsPipeline().push_back(
513  new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
514  pipeline_mng->getOpDomainRhsPipeline().push_back(
515  new OpSetHOInvJacToScalarBases(space, inv_jac_ptr));
516  }
517 
518  pipeline_mng->getOpDomainRhsPipeline().push_back(
519  new OpValsDiffVals(vals, diff_vals, true));
520  pipeline_mng->getOpDomainRhsPipeline().push_back(
521  new OpCalculateScalarFieldValues("FIELD1", ptr_values));
522  pipeline_mng->getOpDomainRhsPipeline().push_back(
524  ptr_diff_vals));
525  pipeline_mng->getOpDomainRhsPipeline().push_back(new OpCheckValsDiffVals(
526  vals, diff_vals, ptr_values, ptr_diff_vals, true));
527 
528  CHKERR pipeline_mng->loopFiniteElements();
529 
531  };
532 
533  CHKERR assemble_matrices_and_vectors();
534  CHKERR solve_problem();
535  CHKERR check_solution();
536  }
537  CATCH_ERRORS;
538 
540 }
static Index< 'o', 3 > o
static const double eps
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
FieldApproximationBase
approximation base
Definition: definitions.h:71
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:75
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:79
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:77
FieldSpace
approximation spaces
Definition: definitions.h:95
@ L2
field with C-1 continuity
Definition: definitions.h:101
@ H1
continuous field
Definition: definitions.h:98
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ MOFEM_IMPOSIBLE_CASE
Definition: definitions.h:48
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:53
#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
auto integration_rule
auto smartCreateDMVector
Get smart vector from DM.
Definition: DMMoFEM.hpp:956
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:481
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:364
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:311
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:342
FTensor::Index< 'i', SPACE_DIM > i
double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
CoreTmp< 0 > Core
Definition: Core.hpp:1096
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:149
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
const double D
diffusivity
const double r
rate factor
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:24
int main(int argc, char *argv[])
EntitiesFieldData::EntData EntData
static int approx_order
static char help[]
constexpr int SPACE_DIM
static FTensor::Tensor1< double, 3 > fUn(const double x, const double y)
static FTensor::Tensor2< double, 3, 2 > diffFun(const double x, const double y)
static double fUn(const double x, const double y, double z)
static FTensor::Tensor1< double, 2 > diffFun(const double x, const double y, double z)
static FTensor::Tensor1< double, 3 > diffFun(const double x, const double y, double z)
static double fUn(const double x, const double y, double z)
PipelineManager::FaceEle DomainEle
Core (interface) class.
Definition: Core.hpp:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FieldData()
Resturn scalar files as a FTensor of rank 0.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:279
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:323
Calculate jacobian on Hex or other volume which is not simplex.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Set inverse jacobian to base functions.
Set inverse jacobian to base functions.
Modify integration weights on face to take in account higher-order geometry.
Transform local reference derivatives of shape functions to global derivatives.
PipelineManager interface.
MoFEM::FaceElementForcesAndSourcesCore FaceEle
Simple interface for fast problem set-up.
Definition: Simple.hpp:33
MoFEMErrorCode addDomainField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
Definition: Simple.cpp:293
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:216
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:715
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:230
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:508
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:668
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
OpCheckValsDiffVals(VectorDouble &vals, MatrixDouble &diff_vals, boost::shared_ptr< VectorDouble > &ptr_vals, boost::shared_ptr< MatrixDouble > &ptr_diff_vals, bool check_grads)
boost::shared_ptr< VectorDouble > ptrVals
FTensor::Index< 'i', SPACE_DIM > i
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< MatrixDouble > ptrDiffVals
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpValsDiffVals(VectorDouble &vals, MatrixDouble &diff_vals, bool check_grads)
FTensor::Index< 'i', SPACE_DIM > i