v0.13.1
approx_sphere.cpp
Go to the documentation of this file.
1/**
2 * \file approx_sphere.cpp
3 * \example approx_sphere.cpp
4 *
5 */
6
7#include <MoFEM.hpp>
8
9using namespace MoFEM;
10
11static char help[] = "...\n\n";
12
13#include <BasicFiniteElements.hpp>
14
15constexpr int FM_DIM = 2;
16
17template <int DIM> struct ElementsAndOps {};
18
19template <> struct ElementsAndOps<2> {
21};
22
26
29
30constexpr double a = 1;
31constexpr double a2 = a * a;
32constexpr double a4 = a2 * a2;
33
34constexpr double A = 6371220;
35
39
40auto res_J = [](const double x, const double y, const double z) {
41 const double res = (x * x + y * y + z * z - a2);
42 return res;
43};
44
45auto res_J_dx = [](const double x, const double y, const double z) {
46 const double res = res_J(x, y, z);
47 return FTensor::Tensor1<double, 3>{res * (2 * x), res * (2 * y),
48 res * (2 * z)};
49};
50
51auto lhs_J_dx2 = [](const double x, const double y, const double z) {
52 const double res = res_J(x, y, z);
54
55 (res * 2 + (4 * x * x)),
56 (4 * y * x),
57 (4 * z * x),
58
59 (4 * x * y),
60 (2 * res + (4 * y * y)),
61 (4 * z * y),
62
63 (4 * x * z),
64 (4 * y * z),
65 (2 * res + (4 * z * z))};
66};
67
68struct OpRhs : public AssemblyDomainEleOp {
69
70 OpRhs(const std::string field_name, boost::shared_ptr<MatrixDouble> x_ptr,
71 boost::shared_ptr<MatrixDouble> dot_x_ptr)
73 xPtr(x_ptr), xDotPtr(dot_x_ptr) {}
74
77
78 auto t_w = getFTensor0IntegrationWeight();
79 auto t_row_base = row_data.getFTensor0N();
80
81 auto t_x0 = getFTensor1CoordsAtGaussPts();
82 auto t_x = getFTensor1FromMat<3>(*xPtr);
83 auto t_dot_x = getFTensor1FromMat<3>(*xDotPtr);
84 auto t_normal = getFTensor1NormalsAtGaussPts();
85
86 constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
87
88 for (int gg = 0; gg != nbIntegrationPts; gg++) {
89
90 constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
91 FTensor::Tensor1<double, 3> t_n{t_x0(0), t_x0(1), t_x0(2)};
92 t_n.normalize();
94 t_P(i, j) = t_n(i) * t_n(j);
95 t_Q(i, j) = t_kd(i, j) - t_P(i, j);
96
97 auto t_J_res = res_J_dx(t_x(0), t_x(1), t_x(2));
98
99 const double alpha = t_w;
100 auto t_nf = getFTensor1FromArray<3, 3>(locF);
101 double l = std::sqrt(t_normal(i) * t_normal(i));
102
104 t_res(i) =
105 alpha * l * ((t_P(i, k) * t_J_res(k) + t_Q(i, k) * t_dot_x(k)));
106
107 int rr = 0;
108 for (; rr != nbRows / 3; ++rr) {
109
110 t_nf(j) += t_row_base * t_res(j);
111
112 ++t_row_base;
113 ++t_nf;
114 }
115 for (; rr < nbRowBaseFunctions; ++rr) {
116 ++t_row_base;
117 }
118
119 ++t_w;
120 ++t_x;
121 ++t_dot_x;
122 ++t_x0;
123 ++t_normal;
124 }
125
127 }
128
129private:
130 boost::shared_ptr<MatrixDouble> xPtr;
131 boost::shared_ptr<MatrixDouble> xDotPtr;
132};
133
134struct OpLhs : public AssemblyDomainEleOp {
135
136 OpLhs(const std::string field_name, boost::shared_ptr<MatrixDouble> x_ptr,
137 boost::shared_ptr<MatrixDouble> dot_x_ptr)
139 AssemblyDomainEleOp::OPROWCOL),
140 xPtr(x_ptr), xDotPtr(dot_x_ptr) {
141 this->sYmm = false;
142 }
143
145 EntitiesFieldData::EntData &col_data) {
147
148 auto t_w = getFTensor0IntegrationWeight();
149 auto t_row_base = row_data.getFTensor0N();
150
151 auto t_x0 = getFTensor1CoordsAtGaussPts();
152 auto t_x = getFTensor1FromMat<3>(*xPtr);
153 auto t_dot_x = getFTensor1FromMat<3>(*xDotPtr);
154 auto t_normal = getFTensor1NormalsAtGaussPts();
155
156 auto get_t_mat = [&](const int rr) {
158 &locMat(rr + 0, 0), &locMat(rr + 0, 1), &locMat(rr + 0, 2),
159
160 &locMat(rr + 1, 0), &locMat(rr + 1, 1), &locMat(rr + 1, 2),
161
162 &locMat(rr + 2, 0), &locMat(rr + 2, 1), &locMat(rr + 2, 2)};
163 };
164
165 constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
166 const double ts_a = getTSa();
167
168 for (int gg = 0; gg != nbIntegrationPts; gg++) {
169
170 constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
171 FTensor::Tensor1<double, 3> t_n{t_x0(0), t_x0(1), t_x0(2)};
172 t_n.normalize();
174 t_P(i, j) = t_n(i) * t_n(j);
175 t_Q(i, j) = t_kd(i, j) - t_P(i, j);
176
177 auto t_J_lhs = lhs_J_dx2(t_x(0), t_x(1), t_x(2));
178 double l = std::sqrt(t_normal(i) * t_normal(i));
179
180 const double alpha = t_w;
182 t_lhs(i, j) =
183 (alpha * l) * (t_P(i, k) * t_J_lhs(k, j) + t_Q(i, j) * ts_a);
184
185 int rr = 0;
186 for (; rr != nbRows / 3; rr++) {
187
188 auto t_col_base = col_data.getFTensor0N(gg, 0);
189 auto t_mat = get_t_mat(3 * rr);
190
191 for (int cc = 0; cc != nbCols / 3; cc++) {
192
193 const double rc = t_row_base * t_col_base;
194 t_mat(i, j) += rc * t_lhs(i, j);
195
196 ++t_col_base;
197 ++t_mat;
198 }
199 ++t_row_base;
200 }
201
202 for (; rr < nbRowBaseFunctions; ++rr)
203 ++t_row_base;
204
205 ++t_w;
206 ++t_x;
207 ++t_x0;
208 ++t_normal;
209 }
210
212 }
213
214private:
215 boost::shared_ptr<MatrixDouble> xPtr;
216 boost::shared_ptr<MatrixDouble> xDotPtr;
217};
218
219struct OpError : public DomainEleOp {
220
221 OpError(const std::string field_name, boost::shared_ptr<MatrixDouble> x_ptr)
223 xPtr(x_ptr) {
224
225 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
226 }
227
229
231
232 auto t_w = getFTensor0IntegrationWeight();
233 auto t_x = getFTensor1FromMat<3>(*xPtr);
234 auto t_normal = getFTensor1NormalsAtGaussPts();
235 auto nb_integration_pts = getGaussPts().size2();
236
237 double error = 0;
238
239 for (int gg = 0; gg != nb_integration_pts; gg++) {
240
241 double l = std::sqrt(t_normal(i) * t_normal(i));
242 error += t_w * l * std::abs((t_x(i) * t_x(i) - A * A));
243
244 ++t_w;
245 ++t_x;
246 ++t_normal;
247 }
248
249 CHKERR VecSetValue(errorVec, 0, error, ADD_VALUES);
250
252 }
253
255
256private:
257 boost::shared_ptr<MatrixDouble> xPtr;
258};
259
261
263
264 ApproxSphere(MoFEM::Interface &m_field) : mField(m_field) {}
265
267
268private:
270
277};
278
279//! [Run programme]
285 CHKERR setOPs();
289}
290//! [Run programme]
291
295}
296
297//! [Read mesh]
301 CHKERR simple->getOptions();
302 CHKERR simple->loadFile();
304}
305//! [Read mesh]
306
307//! [Set up problem]
310
312 CHKERR simple->addDomainField("HO_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE, 3);
313 CHKERR simple->addDataField("HO_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE, 3);
314
315 int order = 3;
316 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
317 CHKERR simple->setFieldOrder("HO_POSITIONS", order);
318 CHKERR simple->setUp();
319
321}
322//! [Set up problem]
323
324//! [Push operators to pipeline]
328
329 auto integration_rule = [](int, int, int approx_order) {
330 return 3 * approx_order;
331 };
334
335 auto x_ptr = boost::make_shared<MatrixDouble>();
336 auto dot_x_ptr = boost::make_shared<MatrixDouble>();
337 auto det_ptr = boost::make_shared<VectorDouble>();
338 auto jac_ptr = boost::make_shared<MatrixDouble>();
339 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
340
341 auto def_ops = [&](auto &pipeline) {
342 pipeline.push_back(
343 new OpCalculateVectorFieldValues<3>("HO_POSITIONS", x_ptr));
344 pipeline.push_back(
345 new OpCalculateVectorFieldValuesDot<3>("HO_POSITIONS", dot_x_ptr));
346 };
347
348 def_ops(pipeline_mng->getOpDomainRhsPipeline());
349 def_ops(pipeline_mng->getOpDomainLhsPipeline());
350
351 pipeline_mng->getOpDomainRhsPipeline().push_back(
352 new OpRhs("HO_POSITIONS", x_ptr, dot_x_ptr));
353 pipeline_mng->getOpDomainLhsPipeline().push_back(
354 new OpLhs("HO_POSITIONS", x_ptr, dot_x_ptr));
355
357}
358//! [Push operators to pipeline]
359
360//! [Solve]
363
364 // Project HO geometry from mesh
365 Projection10NodeCoordsOnField ent_method_material(mField, "HO_POSITIONS");
366 CHKERR mField.loop_dofs("HO_POSITIONS", ent_method_material);
367
368 auto *simple = mField.getInterface<Simple>();
369 auto *pipeline_mng = mField.getInterface<PipelineManager>();
370
371 auto dm = simple->getDM();
373 ts = pipeline_mng->createTSIM();
374
375 double ftime = 1;
376 CHKERR TSSetMaxSteps(ts, 1);
377 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
378
379 auto T = smartCreateDMVector(simple->getDM());
380 CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
381 SCATTER_FORWARD);
382 CHKERR TSSetSolution(ts, T);
383 CHKERR TSSetFromOptions(ts);
384
385 CHKERR TSSolve(ts, NULL);
386 CHKERR TSGetTime(ts, &ftime);
387
388 CHKERR mField.getInterface<FieldBlas>()->fieldScale(A, "HO_POSITIONS");
389
391}
392
393//! [Solve]
396
397 auto x_ptr = boost::make_shared<MatrixDouble>();
398 auto det_ptr = boost::make_shared<VectorDouble>();
399 auto jac_ptr = boost::make_shared<MatrixDouble>();
400 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
401
403 auto dm = simple->getDM();
404
405 auto post_proc_fe =
406 boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(mField);
407
408 post_proc_fe->getOpPtrVector().push_back(
409 new OpCalculateVectorFieldValues<3>("HO_POSITIONS", x_ptr));
410
412
413 post_proc_fe->getOpPtrVector().push_back(
414
415 new OpPPMap(post_proc_fe->getPostProcMesh(),
416 post_proc_fe->getMapGaussPts(),
417
418 {},
419
420 {{"HO_POSITIONS", x_ptr}},
421
422 {}, {}
423
424 )
425
426 );
427
428 CHKERR DMoFEMLoopFiniteElements(dm, "dFE", post_proc_fe);
429 CHKERR post_proc_fe->writeFile("out_approx.h5m");
430
431 auto error_fe = boost::make_shared<DomainEle>(mField);
432
433 error_fe->getOpPtrVector().push_back(
434 new OpGetHONormalsOnFace("HO_POSITIONS"));
435 error_fe->getOpPtrVector().push_back(
436 new OpCalculateVectorFieldValues<3>("HO_POSITIONS", x_ptr));
437 error_fe->getOpPtrVector().push_back(new OpError("HO_POSITIONS", x_ptr));
438
439 error_fe->preProcessHook = [&]() {
441 MOFEM_LOG("EXAMPLE", Sev::inform) << "Create vec ";
443 mField.get_comm(), (!mField.get_comm_rank()) ? 1 : 0, 1);
444 VecZeroEntries(OpError::errorVec);
446 };
447
448 error_fe->postProcessHook = [&]() {
450 CHKERR VecAssemblyBegin(OpError::errorVec);
451 CHKERR VecAssemblyEnd(OpError::errorVec);
452 double error2;
453 CHKERR VecSum(OpError::errorVec, &error2);
454 MOFEM_LOG("EXAMPLE", Sev::inform)
455 << "Error " << std::sqrt(error2 / (4 * M_PI * A * A));
456 OpError::errorVec.reset();
458 };
459
460 CHKERR DMoFEMLoopFiniteElements(dm, "dFE", error_fe);
461
462 CHKERR simple->deleteDM();
463 CHKERR simple->deleteFiniteElements();
464 if (mField.get_comm_size() > 1)
465 CHKERR mField.get_moab().write_file("out_ho_mesh.h5m", "MOAB",
466 "PARALLEL=WRITE_PART");
467 else
468 CHKERR mField.get_moab().write_file("out_ho_mesh.h5m");
470}
471//! [Postprocess results]
472
473int main(int argc, char *argv[]) {
474
475 // Initialisation of MoFEM/PETSc and MOAB data structures
476 const char param_file[] = "param_file.petsc";
478
479 auto core_log = logging::core::get();
480 core_log->add_sink(
481 LogManager::createSink(LogManager::getStrmWorld(), "EXAMPLE"));
482 LogManager::setLog("EXAMPLE");
483 MOFEM_LOG_TAG("EXAMPLE", "example");
484
485 try {
486
487 //! [Register MoFEM discrete manager in PETSc]
488 DMType dm_name = "DMMOFEM";
489 CHKERR DMRegister_MoFEM(dm_name);
490 //! [Register MoFEM discrete manager in PETSc
491
492 //! [Create MoAB]
493 moab::Core mb_instance; ///< mesh database
494 moab::Interface &moab = mb_instance; ///< mesh database interface
495 //! [Create MoAB]
496
497 //! [Create MoFEM]
498 MoFEM::Core core(moab); ///< finite element database
499 MoFEM::Interface &m_field = core; ///< finite element database insterface
500 //! [Create MoFEM]
501
502 //! [ApproxSphere]
503 ApproxSphere ex(m_field);
504 CHKERR ex.runProblem();
505 //! [ApproxSphere]
506 }
508
510}
std::string param_file
ForcesAndSourcesCore::UserDataOperator UserDataOperator
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
constexpr double a2
int main(int argc, char *argv[])
[Postprocess results]
EntitiesFieldData::EntData EntData
constexpr double a
static char help[]
constexpr double a4
constexpr double A
auto lhs_J_dx2
constexpr int FM_DIM
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
auto res_J_dx
auto res_J
FTensor::Index< 'k', 3 > k
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Kronecker Delta class.
Tensor1< T, Tensor_Dim > normalize()
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ 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
#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
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
auto integration_rule
constexpr auto t_kd
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
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:533
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:965
boost::ptr_vector< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
boost::ptr_vector< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:332
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
FTensor::Index< 'l', 3 > l
const double T
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.
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
CoreTmp< 0 > Core
Definition: Core.hpp:1086
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
constexpr auto field_name
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
MoFEMErrorCode readMesh()
[Read mesh]
MoFEMErrorCode solveSystem()
[Push operators to pipeline]
MoFEM::Interface & mField
MoFEMErrorCode setupProblem()
[Read mesh]
ApproxSphere(MoFEM::Interface &m_field)
MoFEMErrorCode runProblem()
[Run programme]
MoFEMErrorCode setOPs()
[Set up problem]
MoFEMErrorCode getOptions()
[Run programme]
MoFEMErrorCode outputResults()
[Solve]
virtual int get_comm_size() const =0
virtual moab::Interface & get_moab()=0
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::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Basic algebra on fields.
Definition: FieldBlas.hpp:21
auto getFTensor0IntegrationWeight()
Get integration weights.
@ OPROW
operator doWork function is executed on FE rows
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
VectorDouble locF
local entity vector
int nbRows
number of dofs on rows
int nbIntegrationPts
number of integration points
MatrixDouble locMat
local entity block matrix
int nbCols
number if dof on column
int nbRowBaseFunctions
number or row base functions
Approximate field valuse for given petsc vector.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Calculate normals at Gauss points of triangle element.
Post post-proc data at points from hash maps.
PipelineManager interface.
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Projection of edge entities with one mid-node on hierarchical basis.
Simple interface for fast problem set-up.
Definition: Simple.hpp:26
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
OpError(const std::string field_name, boost::shared_ptr< MatrixDouble > x_ptr)
static SmartPetscObj< Vec > errorVec
boost::shared_ptr< MatrixDouble > xPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate grad-grad operator.
boost::shared_ptr< MatrixDouble > xPtr
boost::shared_ptr< MatrixDouble > xDotPtr
OpLhs(const std::string field_name, boost::shared_ptr< MatrixDouble > x_ptr, boost::shared_ptr< MatrixDouble > dot_x_ptr)
boost::shared_ptr< MatrixDouble > xPtr
boost::shared_ptr< MatrixDouble > xDotPtr
OpRhs(const std::string field_name, boost::shared_ptr< MatrixDouble > x_ptr, boost::shared_ptr< MatrixDouble > dot_x_ptr)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data)
Class dedicated to integrate operator.