v0.13.1
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 }
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 auto t_row_base = data.getFTensor0N();
150 auto t_diff_row_base = data.getFTensor1DiffN<2>();
151
152 std::array<double, 3> error = {0, 0,
153 0}; // array for storing operator errors
154
155 for (int gg = 0; gg != nb_integration_pts; ++gg) {
156
157 const double alpha = t_w * volume;
158
159 // Calculate errors
160
161 double diff = t_val - fun(t_coords(0), t_coords(1), t_coords(2));
162 error[0] += alpha * pow(diff, 2);
163 auto t_diff_grad = diff_fun(t_coords(0), t_coords(1), t_coords(2));
164 t_diff_grad(i) -= t_grad_val(i);
165
166 error[1] += alpha * t_diff_grad(i) *
167 t_diff_grad(i); // note push forward derivatives
168
170 MOFEM_LOG("SELF", Sev::noisy) << "t_hessian_val " << t_hessian_push;
171
172 // hessian expected to have symmetry
173 if (std::abs(t_hessian_val(0, 1) - t_hessian_val(1, 0)) >
174 std::numeric_limits<float>::epsilon()) {
175 MOFEM_LOG("SELF", Sev::error) << "t_hessian_val " << t_hessian_val;
176 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
177 "Hessian should be symmetric");
178 }
179
180 auto t_diff_hessian = diff2_fun(t_coords(0), t_coords(1), t_coords(2));
181 t_diff_hessian(i, j) -= t_hessian_val(i, j);
182 error[2] = t_diff_hessian(i, j) * t_diff_hessian(i, j);
183
184 // move data to next integration point
185 ++t_w;
186 ++t_val;
187 ++t_grad_val;
188 ++t_hessian_val;
189 ++t_inv_jac;
190 ++t_coords;
191 }
192
193 // assemble error vector
194 std::array<int, 3> index = {0, 1, 2};
195 CHKERR VecSetValues(commonDataPtr->L2Vec, 3, index.data(), error.data(),
196 ADD_VALUES);
197
199 }
200};
201
202//! [Run programme]
211}
212//! [Run programme]
213
214//! [Read mesh]
217
221
223}
224//! [Read mesh]
225
226//! [Set up problem]
229 // Add field
232 constexpr int order = 4;
235
237}
238//! [Set up problem]
239
240//! [Push operators to pipeline]
243
244 auto rule = [](int, int, int p) -> int { return 2 * p; };
245
247 CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule);
248 CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
249
250 auto beta = [](const double, const double, const double) { return 1; };
251 pipeline_mng->getOpDomainLhsPipeline().push_back(
253 pipeline_mng->getOpDomainRhsPipeline().push_back(
255
257}
258//! [Push operators to pipeline]
259
260//! [Solve]
264
265 MOFEM_LOG("WORLD", Sev::inform) << "Solve problem";
266
267 auto solver = pipeline_mng->createKSP();
268 CHKERR KSPSetFromOptions(solver);
269 CHKERR KSPSetUp(solver);
270
271 auto dm = simpleInterface->getDM();
272 auto D = smartCreateDMVector(dm);
273 auto F = smartVectorDuplicate(D);
274
275 CHKERR KSPSolve(solver, F, D);
276 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
277 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
278 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
280}
281
282//! [Check results]
286 pipeline_mng->getDomainLhsFE().reset();
287 pipeline_mng->getDomainRhsFE().reset();
288 pipeline_mng->getOpDomainRhsPipeline().clear();
289
290 auto rule = [](int, int, int p) -> int { return 2 * p; };
292 rule); // set integration rule
293
294 // create data structures for operator
295 auto common_data_ptr = boost::make_shared<CommonData>();
296 common_data_ptr->L2Vec = createSmartVectorMPI(
297 mField.get_comm(), (!mField.get_comm_rank()) ? 3 : 0, 3);
298 common_data_ptr->approxVals = boost::make_shared<VectorDouble>();
299 common_data_ptr->approxGradVals = boost::make_shared<MatrixDouble>();
300 common_data_ptr->approxHessianVals = boost::make_shared<MatrixDouble>();
301
302 // create data strutires for evaluation of higher order derivatives
303 auto base_mass = boost::make_shared<MatrixDouble>();
304 auto data_l2 = boost::make_shared<EntitiesFieldData>(MBENTITYSET);
305 auto jac_ptr = boost::make_shared<MatrixDouble>();
306 auto det_ptr = boost::make_shared<VectorDouble>();
307 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
308 common_data_ptr->invJacPtr = inv_jac_ptr;
309
310 // calculate jacobian at integration points
311 pipeline_mng->getOpDomainRhsPipeline().push_back(
312 new OpCalculateHOJacForFace(jac_ptr));
313 // calculate incerse of jacobian at integration points
314 pipeline_mng->getOpDomainRhsPipeline().push_back(
315 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
316
317 // calculate mass matrix to project derivatives
318 pipeline_mng->getOpDomainRhsPipeline().push_back(
319 new OpBaseDerivativesMass<BASE_DIM>(base_mass, data_l2,
321 // calculate second derivative of base functions, i.e. hessian
322 pipeline_mng->getOpDomainRhsPipeline().push_back(
323 new OpBaseDerivativesNext<BASE_DIM>(BaseDerivatives::SecondDerivative,
324 base_mass, data_l2,
326 // calculate third derivative
327 pipeline_mng->getOpDomainRhsPipeline().push_back(
328 new OpBaseDerivativesNext<BASE_DIM>(BaseDerivatives::ThirdDerivative,
329 base_mass, data_l2,
331 // calculate forth derivative
332 pipeline_mng->getOpDomainRhsPipeline().push_back(
333 new OpBaseDerivativesNext<BASE_DIM>(BaseDerivatives::ForthDerivative,
334 base_mass, data_l2,
336
337 // push first base directives tp physical element shape
338 pipeline_mng->getOpDomainRhsPipeline().push_back(
339 new OpSetInvJacH1ForFace<1>(inv_jac_ptr));
340 // push second base directives tp physical element shape
341 pipeline_mng->getOpDomainRhsPipeline().push_back(
342 new OpSetInvJacH1ForFace<2>(inv_jac_ptr));
343
344 // calculate value of function at integration points
345 pipeline_mng->getOpDomainRhsPipeline().push_back(
347 common_data_ptr->approxVals));
348
349 // calculate gradient at integration points
350 pipeline_mng->getOpDomainRhsPipeline().push_back(
352 FIELD_NAME, common_data_ptr->approxGradVals));
353
354
355
356 // calculate hessian at integration points
357 pipeline_mng->getOpDomainRhsPipeline().push_back(
359 FIELD_NAME, common_data_ptr->approxHessianVals));
360
361 //FIXME: Note third and forth derivative is calculated but not properly tested
362
363 // push operator to evaluate errors
364 pipeline_mng->getOpDomainRhsPipeline().push_back(
365 new OpError(common_data_ptr));
366
367 // zero error vector and iterate over all elements on the mesh
368 CHKERR VecZeroEntries(common_data_ptr->L2Vec);
369 CHKERR pipeline_mng->loopFiniteElements();
370
371 // assemble error vector
372 CHKERR VecAssemblyBegin(common_data_ptr->L2Vec);
373 CHKERR VecAssemblyEnd(common_data_ptr->L2Vec);
374
375 // check if errors are small and print results
376
377 constexpr double eps = 1e-8;
378
379 const double *array;
380 CHKERR VecGetArrayRead(common_data_ptr->L2Vec, &array);
381 MOFEM_LOG_C("WORLD", Sev::inform,
382 "Error %6.4e Diff Error %6.4e Diff2 Error %6.4e\n",
383 std::sqrt(array[0]), std::sqrt(array[1]), std::sqrt(array[2]));
384 if (std::sqrt(array[0]) > eps)
385 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "Wrong function value");
386 if (std::sqrt(array[1]) > eps)
387 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
388 "Wrong function first direcative");
389 if (std::sqrt(array[2]) > eps)
390 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
391 "Wrong function second direcative");
392
393 CHKERR VecRestoreArrayRead(common_data_ptr->L2Vec, &array);
394
396}
397//! [Check results]
398
399int main(int argc, char *argv[]) {
400
401 // Initialisation of MoFEM/PETSc and MOAB data structures
402 MoFEM::Core::Initialize(&argc, &argv, NULL, help);
403
404 try {
405
406 //! [Register MoFEM discrete manager in PETSc]
407 DMType dm_name = "DMMOFEM";
408 CHKERR DMRegister_MoFEM(dm_name);
409 //! [Register MoFEM discrete manager in PETSc
410
411 //! [Create MoAB]
412 moab::Core mb_instance; ///< mesh database
413 moab::Interface &moab = mb_instance; ///< mesh database interface
414 //! [Create MoAB]
415
416 //! [Create MoFEM]
417 MoFEM::Core core(moab); ///< finite element database
418 MoFEM::Interface &m_field = core; ///< finite element database insterface
419 //! [Create MoFEM]
420
421 //! [AtomTest]
422 AtomTest ex(m_field);
423 CHKERR ex.runProblem();
424 //! [AtomTest]
425 }
427
429}
static Index< 'p', 3 > p
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:304
MoFEM::FaceElementForcesAndSourcesCore FaceEle
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
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:470
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:965
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
boost::ptr_vector< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
boost::ptr_vector< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
FTensor::Index< 'i', SPACE_DIM > i
auto diff2_fun
Function second derivative.
int main(int argc, char *argv[])
[Check results]
EntitiesFieldData::EntData EntData
Data on entities.
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.
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: MoFEM.hpp:24
auto createSmartVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
CoreTmp< 0 > Core
Definition: Core.hpp:1086
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
const double D
diffusivity
Collected data use d by operator to evaluate errors for the test.
boost::shared_ptr< VectorDouble > approxVals
boost::shared_ptr< MatrixDouble > approxGradVals
SmartPetscObj< Vec > L2Vec
boost::shared_ptr< MatrixDouble > invJacPtr
boost::shared_ptr< MatrixDouble > approxHessianVals
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)
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.
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()
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Simple interface for fast problem set-up.
Definition: Simple.hpp:26
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:374
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:290
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:806
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:304
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:589
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:749
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.