v0.14.0
Loading...
Searching...
No Matches
higher_derivatives.cpp
Go to the documentation of this file.
1/**
2 * \example higher_derivatives.cpp
3 *
4 * Testing higher derivatives of base functions
5 *
6 */
7
8
9
10#include <MoFEM.hpp>
11
12using namespace MoFEM;
13
14static char help[] = "...\n\n";
15
16constexpr char FIELD_NAME[] = "U";
17constexpr int BASE_DIM = 1;
18constexpr int FIELD_DIM = 1;
19constexpr int SPACE_DIM = 2;
20
21template <int DIM> struct ElementsAndOps {};
22
23template <> struct ElementsAndOps<2> {
26};
27
28using DomainEle = ElementsAndOps<SPACE_DIM>::DomainEle; ///< Finite elenent type
30 DomainEle::UserDataOperator; ///< Finire element operator type
31using EntData = EntitiesFieldData::EntData; ///< Data on entities
32
33/**
34 * @brief Function to approximate
35 *
36 */
37auto fun = [](const double x, const double y, const double z) {
38 return x * x + y * y + x * y + pow(x, 3) + pow(y, 3) + pow(x, 4) + pow(y, 4);
39};
40
41/**
42 * @brief Function derivative
43 *
44 */
45auto diff_fun = [](const double x, const double y, const double z) {
47 2 * x + y + 3 * pow(x, 2) + 4 * pow(x, 3),
48 2 * y + x + 3 * pow(y, 2) + 4 * pow(y, 3)};
49};
50
51/**
52 * @brief Function second derivative
53 *
54 */
55auto diff2_fun = [](const double x, const double y, const double z) {
57 2 + 6 * x + 12 * pow(x, 2), 1.,
58
59 1., 2 + 6 * y + 12 * pow(y, 2)};
60};
61
62/**
63 * @brief OPerator to integrate mass matrix for least square approximation
64 *
65 */
68
69/**
70 * @brief Operator to integrate the right hand side matrix for the problem
71 *
72 */
75
76struct AtomTest {
77
78 AtomTest(MoFEM::Interface &m_field) : mField(m_field) {}
79
81
82private:
85
91
92 /**
93 * @brief Collected data use d by operator to evaluate errors for the test
94 *
95 */
96 struct CommonData {
97 boost::shared_ptr<MatrixDouble> invJacPtr;
98 boost::shared_ptr<VectorDouble> approxVals;
99 boost::shared_ptr<MatrixDouble> approxGradVals;
100 boost::shared_ptr<MatrixDouble> approxHessianVals;
102 };
103
104 /**
105 * @brief Operator to evaluate errors
106 *
107 */
108 struct OpError;
109};
110
111/**
112 * @brief Operator to evaluate errors
113 *
114 */
116 boost::shared_ptr<CommonData> commonDataPtr;
117 OpError(boost::shared_ptr<CommonData> &common_data_ptr)
118 : DomainEleOp(FIELD_NAME, OPROW), commonDataPtr(common_data_ptr) {
119 std::fill(doEntities.begin(), doEntities.end(), false);
120 doEntities[MBVERTEX] = true;
121 }
122 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
124
125 const int nb_integration_pts = getGaussPts().size2();
126 auto t_w = getFTensor0IntegrationWeight(); // ger integration weights
127 auto t_val = getFTensor0FromVec(*(
128 commonDataPtr->approxVals)); // get function approximation at gauss pts
129 auto t_grad_val = getFTensor1FromMat<SPACE_DIM>(
131 ->approxGradVals)); // get gradient of approximation at gauss pts
132 auto t_hessian_val = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(
133 *(commonDataPtr)->approxHessianVals); // get hessian of approximation
134 // at integration pts
135
136 auto t_inv_jac = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(
137 *(commonDataPtr->invJacPtr)); // get inverse of element jacobian
138 auto t_coords = getFTensor1CoordsAtGaussPts(); // get coordinates of
139 // integration points
140
141 // Indices used for tensor operations
146
147 const double volume = getMeasure(); // get finite element area
148
149
150 std::array<double, 3> error = {0, 0,
151 0}; // array for storing operator errors
152
153 for (int gg = 0; gg != nb_integration_pts; ++gg) {
154
155 const double alpha = t_w * volume;
156
157 // Calculate errors
158
159 double diff = t_val - fun(t_coords(0), t_coords(1), t_coords(2));
160 error[0] += alpha * pow(diff, 2);
161 auto t_diff_grad = diff_fun(t_coords(0), t_coords(1), t_coords(2));
162 t_diff_grad(i) -= t_grad_val(i);
163
164 error[1] += alpha * t_diff_grad(i) *
165 t_diff_grad(i); // note push forward derivatives
166
168 MOFEM_LOG("SELF", Sev::noisy) << "t_hessian_val " << t_hessian_push;
169
170 // hessian expected to have symmetry
171 if (std::abs(t_hessian_val(0, 1) - t_hessian_val(1, 0)) >
172 std::numeric_limits<float>::epsilon()) {
173 MOFEM_LOG("SELF", Sev::error) << "t_hessian_val " << t_hessian_val;
174 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
175 "Hessian should be symmetric");
176 }
177
178 auto t_diff_hessian = diff2_fun(t_coords(0), t_coords(1), t_coords(2));
179 t_diff_hessian(i, j) -= t_hessian_val(i, j);
180 error[2] = t_diff_hessian(i, j) * t_diff_hessian(i, j);
181
182 // move data to next integration point
183 ++t_w;
184 ++t_val;
185 ++t_grad_val;
186 ++t_hessian_val;
187 ++t_inv_jac;
188 ++t_coords;
189 }
190
191 // assemble error vector
192 std::array<int, 3> index = {0, 1, 2};
193 CHKERR VecSetValues(commonDataPtr->L2Vec, 3, index.data(), error.data(),
194 ADD_VALUES);
195
197 }
198};
199
200//! [Run programme]
209}
210//! [Run programme]
211
212//! [Read mesh]
215
219
221}
222//! [Read mesh]
223
224//! [Set up problem]
227 // Add field
230 constexpr int order = 4;
233
235}
236//! [Set up problem]
237
238//! [Push operators to pipeline]
241
242 auto rule = [](int, int, int p) -> int { return 2 * p; };
243
245 CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule);
246 CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
247
248 auto beta = [](const double, const double, const double) { return 1; };
249 pipeline_mng->getOpDomainLhsPipeline().push_back(
251 pipeline_mng->getOpDomainRhsPipeline().push_back(
253
255}
256//! [Push operators to pipeline]
257
258//! [Solve]
262
263 MOFEM_LOG("WORLD", Sev::inform) << "Solve problem";
264
265 auto solver = pipeline_mng->createKSP();
266 CHKERR KSPSetFromOptions(solver);
267 CHKERR KSPSetUp(solver);
268
269 auto dm = simpleInterface->getDM();
270 auto D = createDMVector(dm);
271 auto F = vectorDuplicate(D);
272
273 CHKERR KSPSolve(solver, F, D);
274 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
275 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
276 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
278}
279
280//! [Check results]
284 pipeline_mng->getDomainLhsFE().reset();
285 pipeline_mng->getDomainRhsFE().reset();
286 pipeline_mng->getOpDomainRhsPipeline().clear();
287
288 auto rule = [](int, int, int p) -> int { return 2 * p; };
290 rule); // set integration rule
291
292 // create data structures for operator
293 auto common_data_ptr = boost::make_shared<CommonData>();
294 common_data_ptr->L2Vec = createVectorMPI(
295 mField.get_comm(), (!mField.get_comm_rank()) ? 3 : 0, 3);
296 common_data_ptr->approxVals = boost::make_shared<VectorDouble>();
297 common_data_ptr->approxGradVals = boost::make_shared<MatrixDouble>();
298 common_data_ptr->approxHessianVals = boost::make_shared<MatrixDouble>();
299
300 // create data strutires for evaluation of higher order derivatives
301 auto base_mass = boost::make_shared<MatrixDouble>();
302 auto data_l2 = boost::make_shared<EntitiesFieldData>(MBENTITYSET);
303 auto jac_ptr = boost::make_shared<MatrixDouble>();
304 auto det_ptr = boost::make_shared<VectorDouble>();
305 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
306 common_data_ptr->invJacPtr = inv_jac_ptr;
307
308 // calculate jacobian at integration points
309 pipeline_mng->getOpDomainRhsPipeline().push_back(
310 new OpCalculateHOJacForFace(jac_ptr));
311 // calculate incerse of jacobian at integration points
312 pipeline_mng->getOpDomainRhsPipeline().push_back(
313 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
314
315 // calculate mass matrix to project derivatives
316 pipeline_mng->getOpDomainRhsPipeline().push_back(
317 new OpBaseDerivativesMass<BASE_DIM>(base_mass, data_l2,
319 // calculate second derivative of base functions, i.e. hessian
320 pipeline_mng->getOpDomainRhsPipeline().push_back(
321 new OpBaseDerivativesNext<BASE_DIM>(BaseDerivatives::SecondDerivative,
322 base_mass, data_l2,
324 // calculate third derivative
325 pipeline_mng->getOpDomainRhsPipeline().push_back(
326 new OpBaseDerivativesNext<BASE_DIM>(BaseDerivatives::ThirdDerivative,
327 base_mass, data_l2,
329 // calculate forth derivative
330 pipeline_mng->getOpDomainRhsPipeline().push_back(
331 new OpBaseDerivativesNext<BASE_DIM>(BaseDerivatives::ForthDerivative,
332 base_mass, data_l2,
334
335 // push first base directives tp physical element shape
336 pipeline_mng->getOpDomainRhsPipeline().push_back(
337 new OpSetInvJacH1ForFace<1>(inv_jac_ptr));
338 // push second base directives tp physical element shape
339 pipeline_mng->getOpDomainRhsPipeline().push_back(
340 new OpSetInvJacH1ForFace<2>(inv_jac_ptr));
341
342 // calculate value of function at integration points
343 pipeline_mng->getOpDomainRhsPipeline().push_back(
345 common_data_ptr->approxVals));
346
347 // calculate gradient at integration points
348 pipeline_mng->getOpDomainRhsPipeline().push_back(
350 FIELD_NAME, common_data_ptr->approxGradVals));
351
352
353
354 // calculate hessian at integration points
355 pipeline_mng->getOpDomainRhsPipeline().push_back(
357 FIELD_NAME, common_data_ptr->approxHessianVals));
358
359 //FIXME: Note third and forth derivative is calculated but not properly tested
360
361 // push operator to evaluate errors
362 pipeline_mng->getOpDomainRhsPipeline().push_back(
363 new OpError(common_data_ptr));
364
365 // zero error vector and iterate over all elements on the mesh
366 CHKERR VecZeroEntries(common_data_ptr->L2Vec);
367 CHKERR pipeline_mng->loopFiniteElements();
368
369 // assemble error vector
370 CHKERR VecAssemblyBegin(common_data_ptr->L2Vec);
371 CHKERR VecAssemblyEnd(common_data_ptr->L2Vec);
372
373 // check if errors are small and print results
374
375 constexpr double eps = 1e-8;
376
377 const double *array;
378 CHKERR VecGetArrayRead(common_data_ptr->L2Vec, &array);
379 MOFEM_LOG_C("WORLD", Sev::inform,
380 "Error %6.4e Diff Error %6.4e Diff2 Error %6.4e\n",
381 std::sqrt(array[0]), std::sqrt(array[1]), std::sqrt(array[2]));
382 if (std::sqrt(array[0]) > eps)
383 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "Wrong function value");
384 if (std::sqrt(array[1]) > eps)
385 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
386 "Wrong function first direcative");
387 if (std::sqrt(array[2]) > eps)
388 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
389 "Wrong function second direcative");
390
391 CHKERR VecRestoreArrayRead(common_data_ptr->L2Vec, &array);
392
394}
395//! [Check results]
396
397int main(int argc, char *argv[]) {
398
399 // Initialisation of MoFEM/PETSc and MOAB data structures
400 MoFEM::Core::Initialize(&argc, &argv, NULL, help);
401
402 try {
403
404 //! [Register MoFEM discrete manager in PETSc]
405 DMType dm_name = "DMMOFEM";
406 CHKERR DMRegister_MoFEM(dm_name);
407 //! [Register MoFEM discrete manager in PETSc
408
409 //! [Create MoAB]
410 moab::Core mb_instance; ///< mesh database
411 moab::Interface &moab = mb_instance; ///< mesh database interface
412 //! [Create MoAB]
413
414 //! [Create MoFEM]
415 MoFEM::Core core(moab); ///< finite element database
416 MoFEM::Interface &m_field = core; ///< finite element database insterface
417 //! [Create MoFEM]
418
419 //! [AtomTest]
420 AtomTest ex(m_field);
421 CHKERR ex.runProblem();
422 //! [AtomTest]
423 }
425
427}
static Index< 'p', 3 > p
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
int main()
Definition: adol-c_atom.cpp:46
static const double eps
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ H1
continuous field
Definition: definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
constexpr int order
@ F
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:509
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1003
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
boost::ptr_deque< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
boost::ptr_deque< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
FTensor::Index< 'i', SPACE_DIM > i
auto diff2_fun
Function second derivative.
static char help[]
constexpr char FIELD_NAME[]
auto fun
Function to approximate.
constexpr int SPACE_DIM
constexpr int FIELD_DIM
constexpr int BASE_DIM
auto diff_fun
Function derivative.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< BASE_DIM, FIELD_DIM > OpDomainMass
OPerator to integrate mass matrix for least square approximation.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< BASE_DIM, FIELD_DIM > OpDomainSource
Operator to integrate the right hand side matrix for the problem.
double D
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Collected data use d by operator to evaluate errors for the test.
boost::shared_ptr< VectorDouble > approxVals
SmartPetscObj< Vec > L2Vec
boost::shared_ptr< MatrixDouble > invJacPtr
boost::shared_ptr< MatrixDouble > approxHessianVals
boost::shared_ptr< MatrixDouble > approxGradVals
Operator to evaluate errors.
OpError(boost::shared_ptr< CommonData > &common_data_ptr)
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Simple * simpleInterface
MoFEMErrorCode checkResults()
MoFEMErrorCode solveSystem()
AtomTest(MoFEM::Interface &m_field)
MoFEMErrorCode setupProblem()
MoFEMErrorCode readMesh()
MoFEMErrorCode assembleSystem()
MoFEM::Interface & mField
MoFEMErrorCode runProblem()
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
Core (interface) class.
Definition: Core.hpp:82
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
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
auto getFTensor0IntegrationWeight()
Get integration weights.
double getMeasure() const
get measure of element
@ OPROW
operator doWork function is executed on FE rows
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Evaluate field gradient values for scalar field, i.e. gradient is tensor rank 1 (vector),...
Get value at integration points for scalar field.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
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:264
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:667
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:194
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:473
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:611
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.