v0.14.0
Classes | Typedefs | Functions | Variables
operators_tests.cpp File Reference
#include <MoFEM.hpp>

Go to the source code of this file.

Classes

struct  ElementsAndOps< DIM >
 
struct  ElementsAndOps< 2 >
 
struct  ElementsAndOps< 3 >
 

Typedefs

using EntData = EntitiesFieldData::EntData
 
using DomainEle = ElementsAndOps< SPACE_DIM >::DomainEle
 
using DomainEleOp = DomainEle::UserDataOperator
 
using PostProcEle = PostProcBrokenMeshInMoab< DomainEle >
 
template<int FIELD_DIM>
using OpGradGrad = FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpGradGrad< 1, FIELD_DIM, SPACE_DIM >
 
template<int FIELD_DIM>
using OpGradTimesTensor = FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM >
 
template<int FIELD_DIM>
using OpConvectiveTermRhs = FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpConvectiveTermRhs< 1, FIELD_DIM, SPACE_DIM >
 
template<int FIELD_DIM>
using OpConvectiveTermLhsDu = FormsIntegrators< DomainEleOp >::Assembly< A >::TriLinearForm< I >::OpConvectiveTermLhsDu< 1, FIELD_DIM, SPACE_DIM >
 
template<int FIELD_DIM>
using OpConvectiveTermLhsDy = FormsIntegrators< DomainEleOp >::Assembly< A >::TriLinearForm< I >::OpConvectiveTermLhsDy< 1, FIELD_DIM, SPACE_DIM >
 

Functions

int main (int argc, char *argv[])
 

Variables

static char help [] = "...\n\n"
 
constexpr int SPACE_DIM = 3
 [Define dimension] More...
 
constexpr AssemblyType A = AssemblyType::PETSC
 
constexpr IntegrationType I
 
constexpr bool debug = false
 

Typedef Documentation

◆ DomainEle

Definition at line 35 of file operators_tests.cpp.

◆ DomainEleOp

Definition at line 36 of file operators_tests.cpp.

◆ EntData

Examples
operators_tests.cpp.

Definition at line 34 of file operators_tests.cpp.

◆ OpConvectiveTermLhsDu

template<int FIELD_DIM>
using OpConvectiveTermLhsDu = FormsIntegrators<DomainEleOp>::Assembly< A>::TriLinearForm<I>::OpConvectiveTermLhsDu<1, FIELD_DIM, SPACE_DIM>
Examples
operators_tests.cpp, and shallow_wave.cpp.

Definition at line 56 of file operators_tests.cpp.

◆ OpConvectiveTermLhsDy

template<int FIELD_DIM>
using OpConvectiveTermLhsDy = FormsIntegrators<DomainEleOp>::Assembly< A>::TriLinearForm<I>::OpConvectiveTermLhsDy<1, FIELD_DIM, SPACE_DIM>
Examples
operators_tests.cpp, and shallow_wave.cpp.

Definition at line 60 of file operators_tests.cpp.

◆ OpConvectiveTermRhs

template<int FIELD_DIM>
using OpConvectiveTermRhs = FormsIntegrators<DomainEleOp>::Assembly< A>::LinearForm<I>::OpConvectiveTermRhs<1, FIELD_DIM, SPACE_DIM>
Examples
operators_tests.cpp, and shallow_wave.cpp.

Definition at line 52 of file operators_tests.cpp.

◆ OpGradGrad

template<int FIELD_DIM>
using OpGradGrad = FormsIntegrators<DomainEleOp>::Assembly<A>::BiLinearForm< I>::OpGradGrad<1, FIELD_DIM, SPACE_DIM>

Definition at line 44 of file operators_tests.cpp.

◆ OpGradTimesTensor

template<int FIELD_DIM>
using OpGradTimesTensor = FormsIntegrators<DomainEleOp>::Assembly< A>::LinearForm<I>::OpGradTimesTensor<1, FIELD_DIM, SPACE_DIM>

◆ PostProcEle

Definition at line 37 of file operators_tests.cpp.

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)
Examples
operators_tests.cpp.

Definition at line 64 of file operators_tests.cpp.

64  {
65 
66  // initialize petsc
67  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
68 
69  try {
70 
71  // Create MoAB database
72  moab::Core moab_core;
73  moab::Interface &moab = moab_core;
74 
75  // Create MoFEM database and link it to MoAB
76  MoFEM::Core mofem_core(moab);
77  MoFEM::Interface &m_field = mofem_core;
78 
79  // Register DM Manager
80  DMType dm_name = "DMMOFEM";
81  CHKERR DMRegister_MoFEM(dm_name);
82 
83  // Add logging channel for example
84  auto core_log = logging::core::get();
85  core_log->add_sink(
86  LogManager::createSink(LogManager::getStrmWorld(), "OpTester"));
87  LogManager::setLog("OpTester");
88  MOFEM_LOG_TAG("OpTester", "OpTester");
89 
90  // Simple interface
91  auto simple = m_field.getInterface<Simple>();
92 
93  // get options from command line
94  CHKERR simple->getOptions();
95  // load mesh file
96  CHKERR simple->loadFile();
97 
98  // Scalar fields and vector field is tested. Add more fields, i.e. vector
99  // field if needed.
100  CHKERR simple->addDomainField("SCALAR", H1, AINSWORTH_LEGENDRE_BASE, 1);
101  CHKERR simple->addDomainField("VECTOR", H1, AINSWORTH_LEGENDRE_BASE,
102  SPACE_DIM);
103 
104  // set fields order, i.e. for most first cases order is sufficient.
105  CHKERR simple->setFieldOrder("SCALAR", 1);
106  CHKERR simple->setFieldOrder("VECTOR", 1);
107 
108  // setup problem
109  CHKERR simple->setUp();
110 
111  // get operators tester
112  auto opt = m_field.getInterface<OperatorsTester>(); // get interface to
113  // OperatorsTester
114  auto pip = m_field.getInterface<PipelineManager>(); // get interface to
115  // pipeline manager
116 
117  auto post_proc = [&](auto dm, auto f_res, auto out_name) {
119  auto post_proc_fe =
120  boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(m_field);
121 
123 
124  auto scal_vec = boost::make_shared<VectorDouble>();
125  auto vec_mat = boost::make_shared<MatrixDouble>();
126 
127  post_proc_fe->getOpPtrVector().push_back(
128  new OpCalculateScalarFieldValues("SCALAR", scal_vec, f_res));
129  post_proc_fe->getOpPtrVector().push_back(
130  new OpCalculateVectorFieldValues<SPACE_DIM>("VECTOR", vec_mat,
131  f_res));
132 
133  post_proc_fe->getOpPtrVector().push_back(
134 
135  new OpPPMap(
136 
137  post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
138 
139  {{"SCALAR", scal_vec}},
140 
141  {{"VECTOR", vec_mat}},
142 
143  {}, {})
144 
145  );
146 
147  CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
148  post_proc_fe);
149  post_proc_fe->writeFile(out_name);
151  };
152 
153  // Note: I do not push OPs transferring base to physical coordinates. That
154  // is not needed to test OPs, and not required, since we like to test only
155  // operators.
156 
157  // Test grad-grad operator for scalar and vector field specialisation
158  auto TestOpGradGrad = [&]() {
160  MOFEM_LOG("OpTester", Sev::verbose) << "TestOpGradGrad";
161 
162  pip->getOpDomainLhsPipeline().clear();
163  pip->getOpDomainRhsPipeline().clear();
164 
165  pip->setDomainLhsIntegrationRule([](int, int, int o) { return 2 * o; });
166  pip->setDomainLhsIntegrationRule([](int, int, int o) { return 2 * o; });
167 
168  pip->getOpDomainLhsPipeline().push_back(new OpGradGrad<1>(
169  "SCALAR", "SCALAR", [](double, double, double) { return 1; }));
170  pip->getOpDomainLhsPipeline().push_back(new OpGradGrad<SPACE_DIM>(
171  "VECTOR", "VECTOR", [](double, double, double) { return 1; }));
172 
173  auto scl_mat = boost::make_shared<MatrixDouble>();
174  auto vec_mat = boost::make_shared<MatrixDouble>();
175 
176  pip->getOpDomainRhsPipeline().push_back(
177  new OpCalculateScalarFieldGradient<SPACE_DIM>("SCALAR", scl_mat));
178  pip->getOpDomainRhsPipeline().push_back(
180  vec_mat));
181 
182  pip->getOpDomainRhsPipeline().push_back(new OpGradTimesTensor<1>(
183  "SCALAR", scl_mat, [](double, double, double) { return 1; }));
184  pip->getOpDomainRhsPipeline().push_back(new OpGradTimesTensor<SPACE_DIM>(
185  "VECTOR", vec_mat, [](double, double, double) { return 1; }));
186 
187  constexpr double eps = 1e-6;
188 
189  auto x = opt->setRandomFields(simple->getDM(),
190  {{"SCALAR", {-1, 1}}, {"VECTOR", {-1, 1}}});
191  auto diff_x = opt->setRandomFields(
192  simple->getDM(), {{"SCALAR", {-1, 1}}, {"VECTOR", {-1, 1}}});
193 
194  auto diff_res = opt->checkCentralFiniteDifference(
195  simple->getDM(), simple->getDomainFEName(), pip->getDomainRhsFE(),
196  pip->getDomainLhsFE(), x, SmartPetscObj<Vec>(), SmartPetscObj<Vec>(),
197  diff_x, 0, 1, eps);
198 
199  // Calculate norm of difference between directive calculated from finite
200  // difference, and tangent matrix.
201  double fnorm;
202  CHKERR VecNorm(diff_res, NORM_2, &fnorm);
203  MOFEM_LOG_C("OpTester", Sev::inform, "TestOpGradGrad %3.4e", fnorm);
204 
205  constexpr double err = 1e-9;
206  if (fnorm > err)
207  SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
208  "Norm of directional derivative too large");
209 
210  if (debug) {
211  // Example how to plot direction in direction diff_x. If instead
212  // directionalCentralFiniteDifference(...) diff_res is used, then error
213  // on directive is plotted.
214  CHKERR post_proc(simple->getDM(),
215  //
216  opt->directionalCentralFiniteDifference(
217  simple->getDM(), simple->getDomainFEName(),
218  pip->getDomainRhsFE(), x, SmartPetscObj<Vec>(),
219  SmartPetscObj<Vec>(), diff_res, 0, 1, eps),
220  //
221  "out_directional_directive_gradgrad.h5m");
222  }
223 
225  };
226 
227  // test operator for convection (part of material time derivatives)
228  auto TestOpConvection = [&]() {
230  MOFEM_LOG("OpTester", Sev::verbose) << "TestOpConvection";
231 
232  pip->getOpDomainLhsPipeline().clear();
233  pip->getOpDomainRhsPipeline().clear();
234 
235  pip->setDomainLhsIntegrationRule([](int, int, int o) { return 2 * o; });
236  pip->setDomainLhsIntegrationRule([](int, int, int o) { return 2 * o; });
237 
238  auto evaluate_field = [&](auto &pipe) {
239  auto scl_mat = boost::make_shared<MatrixDouble>();
240  auto vec_mat = boost::make_shared<MatrixDouble>();
241  auto dot_vec_vel = boost::make_shared<MatrixDouble>();
242  pipe.push_back(
243  new OpCalculateScalarFieldGradient<SPACE_DIM>("SCALAR", scl_mat));
245  "VECTOR", vec_mat));
247  "VECTOR", dot_vec_vel));
248  return std::make_tuple(scl_mat, vec_mat, dot_vec_vel);
249  };
250 
251  auto [rhs_scl_mat, rhs_vec_mat, rhs_dot_vec_vel] =
252  evaluate_field(pip->getOpDomainRhsPipeline());
253  pip->getOpDomainRhsPipeline().push_back(
254  new OpConvectiveTermRhs<1>("SCALAR", rhs_dot_vec_vel, rhs_scl_mat));
255  pip->getOpDomainRhsPipeline().push_back(
256  new OpConvectiveTermRhs<SPACE_DIM>("VECTOR", rhs_dot_vec_vel,
257  rhs_vec_mat));
258 
259  auto [lhs_scl_mat, lhs_vec_mat, lhs_dot_vec_vel] =
260  evaluate_field(pip->getOpDomainLhsPipeline());
261 
262  // I create op, set scaling function to calculate time directive, and add
263  // operator pinter to pipeline
264  auto op_convective_term_lhs_du_scalar =
265  new OpConvectiveTermLhsDu<1>("SCALAR", "VECTOR", lhs_scl_mat);
266  op_convective_term_lhs_du_scalar->feScalingFun =
267  [](const FEMethod *fe_ptr) { return fe_ptr->ts_a; };
268  pip->getOpDomainLhsPipeline().push_back(op_convective_term_lhs_du_scalar);
269  pip->getOpDomainLhsPipeline().push_back(
270  new OpConvectiveTermLhsDy<1>("SCALAR", "SCALAR", lhs_dot_vec_vel));
271 
272  auto op_convective_term_lhs_du_vec =
273  new OpConvectiveTermLhsDu<SPACE_DIM>("VECTOR", "VECTOR", lhs_vec_mat);
274  op_convective_term_lhs_du_vec->feScalingFun = [](const FEMethod *fe_ptr) {
275  return fe_ptr->ts_a;
276  };
277  pip->getOpDomainLhsPipeline().push_back(op_convective_term_lhs_du_vec);
278  pip->getOpDomainLhsPipeline().push_back(
279  new OpConvectiveTermLhsDy<SPACE_DIM>("VECTOR", "VECTOR",
280  lhs_dot_vec_vel));
281 
282  constexpr double eps = 1e-6;
283 
284  auto x = opt->setRandomFields(simple->getDM(),
285  {{"SCALAR", {-1, 1}}, {"VECTOR", {-1, 1}}});
286  auto x_t = opt->setRandomFields(
287  simple->getDM(), {{"SCALAR", {-1, 1}}, {"VECTOR", {-1, 1}}});
288  auto diff_x = opt->setRandomFields(
289  simple->getDM(), {{"SCALAR", {-1, 1}}, {"VECTOR", {-1, 1}}});
290 
291  auto diff_res = opt->checkCentralFiniteDifference(
292  simple->getDM(), simple->getDomainFEName(), pip->getDomainRhsFE(),
293  pip->getDomainLhsFE(), x, x_t, SmartPetscObj<Vec>(), diff_x, 0, 1,
294  eps);
295 
296  // Calculate norm of difference between directive calculated from finite
297  // difference, and tangent matrix.
298  double fnorm;
299  CHKERR VecNorm(diff_res, NORM_2, &fnorm);
300  MOFEM_LOG_C("OpTester", Sev::inform, "TestOpConvection %3.4e", fnorm);
301 
302  constexpr double err = 1e-9;
303  if (fnorm > err) {
304  SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
305  "Norm of directional derivative too large");
306  }
307 
308  if (debug) {
309  // Example how to plot direction in direction diff_x. If instead
310  // directionalCentralFiniteDifference(...) diff_res is used, then error
311  // on directive is plotted.
312  CHKERR post_proc(simple->getDM(),
313  //
314  opt->directionalCentralFiniteDifference(
315  simple->getDM(), simple->getDomainFEName(),
316  pip->getDomainRhsFE(), x, x_t,
317  SmartPetscObj<Vec>(), diff_res, 0, 1, eps),
318  //
319  "out_directional_directive_convection.h5m");
320  }
321 
323  };
324 
325  CHKERR TestOpGradGrad();
326  CHKERR TestOpConvection();
327  }
328  CATCH_ERRORS;
329 
330  // finish work cleaning memory, getting statistics, etc.
332 
333  return 0;
334 }

Variable Documentation

◆ A

constexpr AssemblyType A = AssemblyType::PETSC
constexpr

◆ debug

constexpr bool debug = false
constexpr
Examples
operators_tests.cpp.

Definition at line 62 of file operators_tests.cpp.

◆ help

char help[] = "...\n\n"
static
Examples
operators_tests.cpp.

Definition at line 16 of file operators_tests.cpp.

◆ I

constexpr IntegrationType I
constexpr

◆ SPACE_DIM

constexpr int SPACE_DIM = 3
constexpr

[Define dimension]

Examples
operators_tests.cpp.

Definition at line 29 of file operators_tests.cpp.

MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:369
OpGradGrad
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpGradGrad< 1, FIELD_DIM, SPACE_DIM > OpGradGrad
Definition: operators_tests.cpp:44
OpConvectiveTermRhs
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpConvectiveTermRhs< 1, FIELD_DIM, SPACE_DIM > OpConvectiveTermRhs
Definition: operators_tests.cpp:52
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:466
OpConvectiveTermLhsDu
FormsIntegrators< DomainEleOp >::Assembly< A >::TriLinearForm< I >::OpConvectiveTermLhsDu< 1, FIELD_DIM, SPACE_DIM > OpConvectiveTermLhsDu
Definition: operators_tests.cpp:56
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1293
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MoFEM::OpCalculateVectorFieldValuesFromPetscVecImpl
Approximate field values for given petsc vector.
Definition: UserDataOperators.hpp:595
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
OpConvectiveTermLhsDy
FormsIntegrators< DomainEleOp >::Assembly< A >::TriLinearForm< I >::OpConvectiveTermLhsDy< 1, FIELD_DIM, SPACE_DIM > OpConvectiveTermLhsDy
Definition: operators_tests.cpp:60
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
OpGradTimesTensor
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
Definition: operators_tests.cpp:48
MoFEM::OperatorsTester
Calculate directional derivative of the right hand side and compare it with tangent matrix derivative...
Definition: OperatorsTester.hpp:21
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1535
debug
constexpr bool debug
Definition: operators_tests.cpp:62
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
MoFEM::SmartPetscObj< Vec >
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:586
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
SPACE_DIM
constexpr int SPACE_DIM
[Define dimension]
Definition: operators_tests.cpp:29
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698
help
static char help[]
Definition: operators_tests.cpp:16