v0.13.0
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 /* This file is part of MoFEM.
8  * MoFEM is free software: you can redistribute it and/or modify it under
9  * the terms of the GNU Lesser General Public License as published by the
10  * Free Software Foundation, either version 3 of the License, or (at your
11  * option) any later version.
12  *
13  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
14  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  * License for more details.
17  *
18  * You should have received a copy of the GNU Lesser General Public
19  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
20 
21 #include <MoFEM.hpp>
22 
23 using namespace MoFEM;
24 
25 static char help[] = "...\n\n";
26 
27 #include <BasicFiniteElements.hpp>
28 
29 constexpr int FM_DIM = 2;
30 
31 template <int DIM> struct ElementsAndOps {};
32 
33 template <> struct ElementsAndOps<2> {
37 };
38 
42 
45 
46 constexpr double a = 1;
47 constexpr double a2 = a * a;
48 constexpr double a4 = a2 * a2;
49 
50 constexpr double A = 6371220;
51 
55 
56 auto res_J = [](const double x, const double y, const double z) {
57  const double res = (x * x + y * y + z * z - a2);
58  return res;
59 };
60 
61 auto res_J_dx = [](const double x, const double y, const double z) {
62  const double res = res_J(x, y, z);
63  return FTensor::Tensor1<double, 3>{res * (2 * x), res * (2 * y),
64  res * (2 * z)};
65 };
66 
67 auto lhs_J_dx2 = [](const double x, const double y, const double z) {
68  const double res = res_J(x, y, z);
70 
71  (res * 2 + (4 * x * x)),
72  (4 * y * x),
73  (4 * z * x),
74 
75  (4 * x * y),
76  (2 * res + (4 * y * y)),
77  (4 * z * y),
78 
79  (4 * x * z),
80  (4 * y * z),
81  (2 * res + (4 * z * z))};
82 };
83 
84 struct OpRhs : public AssemblyDomainEleOp {
85 
86  OpRhs(const std::string field_name, boost::shared_ptr<MatrixDouble> x_ptr,
87  boost::shared_ptr<MatrixDouble> dot_x_ptr)
88  : AssemblyDomainEleOp(field_name, field_name, AssemblyDomainEleOp::OPROW),
89  xPtr(x_ptr), xDotPtr(dot_x_ptr) {}
90 
93 
94  auto t_w = getFTensor0IntegrationWeight();
95  auto t_row_base = row_data.getFTensor0N();
96 
97  auto t_x0 = getFTensor1CoordsAtGaussPts();
98  auto t_x = getFTensor1FromMat<3>(*xPtr);
99  auto t_dot_x = getFTensor1FromMat<3>(*xDotPtr);
100  auto t_normal = getFTensor1NormalsAtGaussPts();
101 
102  constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
103 
104  for (int gg = 0; gg != nbIntegrationPts; gg++) {
105 
106  constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
107  FTensor::Tensor1<double, 3> t_n{t_x0(0), t_x0(1), t_x0(2)};
108  t_n.normalize();
110  t_P(i, j) = t_n(i) * t_n(j);
111  t_Q(i, j) = t_kd(i, j) - t_P(i, j);
112 
113  auto t_J_res = res_J_dx(t_x(0), t_x(1), t_x(2));
114 
115  const double alpha = t_w;
116  auto t_nf = getFTensor1FromArray<3, 3>(locF);
117  double l = std::sqrt(t_normal(i) * t_normal(i));
118 
120  t_res(i) =
121  alpha * l * ((t_P(i, k) * t_J_res(k) + t_Q(i, k) * t_dot_x(k)));
122 
123  int rr = 0;
124  for (; rr != nbRows / 3; ++rr) {
125 
126  t_nf(j) += t_row_base * t_res(j);
127 
128  ++t_row_base;
129  ++t_nf;
130  }
131  for (; rr < nbRowBaseFunctions; ++rr) {
132  ++t_row_base;
133  }
134 
135  ++t_w;
136  ++t_x;
137  ++t_dot_x;
138  ++t_x0;
139  ++t_normal;
140  }
141 
143  }
144 
145 private:
146  boost::shared_ptr<MatrixDouble> xPtr;
147  boost::shared_ptr<MatrixDouble> xDotPtr;
148 };
149 
150 struct OpLhs : public AssemblyDomainEleOp {
151 
152  OpLhs(const std::string field_name, boost::shared_ptr<MatrixDouble> x_ptr,
153  boost::shared_ptr<MatrixDouble> dot_x_ptr)
154  : AssemblyDomainEleOp(field_name, field_name,
155  AssemblyDomainEleOp::OPROWCOL),
156  xPtr(x_ptr), xDotPtr(dot_x_ptr) {
157  this->sYmm = false;
158  }
159 
161  EntitiesFieldData::EntData &col_data) {
163 
164  auto t_w = getFTensor0IntegrationWeight();
165  auto t_row_base = row_data.getFTensor0N();
166 
167  auto t_x0 = getFTensor1CoordsAtGaussPts();
168  auto t_x = getFTensor1FromMat<3>(*xPtr);
169  auto t_dot_x = getFTensor1FromMat<3>(*xDotPtr);
170  auto t_normal = getFTensor1NormalsAtGaussPts();
171 
172  auto get_t_mat = [&](const int rr) {
174  &locMat(rr + 0, 0), &locMat(rr + 0, 1), &locMat(rr + 0, 2),
175 
176  &locMat(rr + 1, 0), &locMat(rr + 1, 1), &locMat(rr + 1, 2),
177 
178  &locMat(rr + 2, 0), &locMat(rr + 2, 1), &locMat(rr + 2, 2)};
179  };
180 
181  constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
182  const double ts_a = getTSa();
183 
184  for (int gg = 0; gg != nbIntegrationPts; gg++) {
185 
186  constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
187  FTensor::Tensor1<double, 3> t_n{t_x0(0), t_x0(1), t_x0(2)};
188  t_n.normalize();
190  t_P(i, j) = t_n(i) * t_n(j);
191  t_Q(i, j) = t_kd(i, j) - t_P(i, j);
192 
193  auto t_J_lhs = lhs_J_dx2(t_x(0), t_x(1), t_x(2));
194  double l = std::sqrt(t_normal(i) * t_normal(i));
195 
196  const double alpha = t_w;
198  t_lhs(i, j) =
199  (alpha * l) * (t_P(i, k) * t_J_lhs(k, j) + t_Q(i, j) * ts_a);
200 
201  int rr = 0;
202  for (; rr != nbRows / 3; rr++) {
203 
204  auto t_col_base = col_data.getFTensor0N(gg, 0);
205  auto t_mat = get_t_mat(3 * rr);
206 
207  for (int cc = 0; cc != nbCols / 3; cc++) {
208 
209  const double rc = t_row_base * t_col_base;
210  t_mat(i, j) += rc * t_lhs(i, j);
211 
212  ++t_col_base;
213  ++t_mat;
214  }
215  ++t_row_base;
216  }
217 
218  for (; rr < nbRowBaseFunctions; ++rr)
219  ++t_row_base;
220 
221  ++t_w;
222  ++t_x;
223  ++t_x0;
224  ++t_normal;
225  }
226 
228  }
229 
230 private:
231  boost::shared_ptr<MatrixDouble> xPtr;
232  boost::shared_ptr<MatrixDouble> xDotPtr;
233 };
234 
235 struct OpError : public DomainEleOp {
236 
237  OpError(const std::string field_name, boost::shared_ptr<MatrixDouble> x_ptr)
238  : DomainEleOp(field_name, field_name, AssemblyDomainEleOp::OPROW),
239  xPtr(x_ptr) {
240 
241  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
242  }
243 
244  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
245 
247 
248  auto t_w = getFTensor0IntegrationWeight();
249  auto t_x = getFTensor1FromMat<3>(*xPtr);
250  auto t_normal = getFTensor1NormalsAtGaussPts();
251  auto nb_integration_pts = getGaussPts().size2();
252 
253  double error = 0;
254 
255  for (int gg = 0; gg != nb_integration_pts; gg++) {
256 
257  double l = std::sqrt(t_normal(i) * t_normal(i));
258  error += t_w * l * std::abs((t_x(i) * t_x(i) - A * A));
259 
260  ++t_w;
261  ++t_x;
262  ++t_normal;
263  }
264 
265  CHKERR VecSetValue(errorVec, 0, error, ADD_VALUES);
266 
268  }
269 
271 
272 private:
273  boost::shared_ptr<MatrixDouble> xPtr;
274 };
275 
277 
278 struct ApproxSphere {
279 
280  ApproxSphere(MoFEM::Interface &m_field) : mField(m_field) {}
281 
282  MoFEMErrorCode runProblem();
283 
284 private:
286 
287  MoFEMErrorCode getOptions();
288  MoFEMErrorCode readMesh();
289  MoFEMErrorCode setupProblem();
290  MoFEMErrorCode setOPs();
291  MoFEMErrorCode solveSystem();
292  MoFEMErrorCode outputResults();
293 };
294 
295 //! [Run programme]
298  CHKERR getOptions();
299  CHKERR readMesh();
300  CHKERR setupProblem();
301  CHKERR setOPs();
302  CHKERR solveSystem();
303  CHKERR outputResults();
305 }
306 //! [Run programme]
307 
311 }
312 
313 //! [Read mesh]
316  auto simple = mField.getInterface<Simple>();
317  CHKERR simple->getOptions();
318  CHKERR simple->loadFile();
320 }
321 //! [Read mesh]
322 
323 //! [Set up problem]
326 
327  auto simple = mField.getInterface<Simple>();
328  CHKERR simple->addDomainField("HO_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE, 3);
329  CHKERR simple->addDataField("HO_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE, 3);
330 
331  int order = 3;
332  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
333  CHKERR simple->setFieldOrder("HO_POSITIONS", order);
334  CHKERR simple->setUp();
335 
337 }
338 //! [Set up problem]
339 
340 //! [Push operators to pipeline]
343  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
344 
345  auto integration_rule = [](int, int, int approx_order) {
346  return 3 * approx_order;
347  };
350 
351  auto x_ptr = boost::make_shared<MatrixDouble>();
352  auto dot_x_ptr = boost::make_shared<MatrixDouble>();
353  auto det_ptr = boost::make_shared<VectorDouble>();
354  auto jac_ptr = boost::make_shared<MatrixDouble>();
355  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
356 
357  auto def_ops = [&](auto &pipeline) {
358  pipeline.push_back(
359  new OpCalculateVectorFieldValues<3>("HO_POSITIONS", x_ptr));
360  pipeline.push_back(
361  new OpCalculateVectorFieldValuesDot<3>("HO_POSITIONS", dot_x_ptr));
362  };
363 
364  def_ops(pipeline_mng->getOpDomainRhsPipeline());
365  def_ops(pipeline_mng->getOpDomainLhsPipeline());
366 
367  pipeline_mng->getOpDomainRhsPipeline().push_back(
368  new OpRhs("HO_POSITIONS", x_ptr, dot_x_ptr));
369  pipeline_mng->getOpDomainLhsPipeline().push_back(
370  new OpLhs("HO_POSITIONS", x_ptr, dot_x_ptr));
371 
373 }
374 //! [Push operators to pipeline]
375 
376 //! [Solve]
379 
380  // Project HO geometry from mesh
381  Projection10NodeCoordsOnField ent_method_material(mField, "HO_POSITIONS");
382  CHKERR mField.loop_dofs("HO_POSITIONS", ent_method_material);
383 
384  auto *simple = mField.getInterface<Simple>();
385  auto *pipeline_mng = mField.getInterface<PipelineManager>();
386 
387  auto dm = simple->getDM();
389  ts = pipeline_mng->createTSIM();
390 
391  double ftime = 1;
392  CHKERR TSSetMaxSteps(ts, 1);
393  CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
394 
395  auto T = smartCreateDMVector(simple->getDM());
396  CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
397  SCATTER_FORWARD);
398  CHKERR TSSetSolution(ts, T);
399  CHKERR TSSetFromOptions(ts);
400 
401  CHKERR TSSolve(ts, NULL);
402  CHKERR TSGetTime(ts, &ftime);
403 
404  CHKERR mField.getInterface<FieldBlas>()->fieldScale(A, "HO_POSITIONS");
405 
407 }
408 
409 //! [Solve]
412 
413  auto simple = mField.getInterface<Simple>();
414  auto dm = simple->getDM();
415 
416  auto post_proc_fe = boost::make_shared<PostProcFaceOnRefinedMesh>(mField);
417  post_proc_fe->generateReferenceElementMesh();
418  post_proc_fe->addFieldValuesPostProc("HO_POSITIONS");
419  CHKERR DMoFEMLoopFiniteElements(dm, "dFE", post_proc_fe);
420  CHKERR post_proc_fe->writeFile("out_approx.h5m");
421 
422  auto error_fe = boost::make_shared<DomainEle>(mField);
423 
424  auto x_ptr = boost::make_shared<MatrixDouble>();
425  auto det_ptr = boost::make_shared<VectorDouble>();
426  auto jac_ptr = boost::make_shared<MatrixDouble>();
427  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
428 
429  error_fe->getOpPtrVector().push_back(
430  new OpGetHONormalsOnFace("HO_POSITIONS"));
431  error_fe->getOpPtrVector().push_back(
432  new OpCalculateVectorFieldValues<3>("HO_POSITIONS", x_ptr));
433  error_fe->getOpPtrVector().push_back(new OpError("HO_POSITIONS", x_ptr));
434 
435  error_fe->preProcessHook = [&]() {
437  MOFEM_LOG("EXAMPLE", Sev::inform) << "Create vec ";
439  mField.get_comm(), (!mField.get_comm_rank()) ? 1 : 0, 1);
440  VecZeroEntries(OpError::errorVec);
442  };
443 
444  error_fe->postProcessHook = [&]() {
446  CHKERR VecAssemblyBegin(OpError::errorVec);
447  CHKERR VecAssemblyEnd(OpError::errorVec);
448  double error2;
449  CHKERR VecSum(OpError::errorVec, &error2);
450  MOFEM_LOG("EXAMPLE", Sev::inform)
451  << "Error " << std::sqrt(error2 / (4 * M_PI * A * A));
452  OpError::errorVec.reset();
454  };
455 
456  CHKERR DMoFEMLoopFiniteElements(dm, "dFE", error_fe);
457 
458  CHKERR simple->deleteDM();
459  CHKERR simple->deleteFiniteElements();
460  if (mField.get_comm_size() > 1)
461  CHKERR mField.get_moab().write_file("out_ho_mesh.h5m", "MOAB",
462  "PARALLEL=WRITE_PART");
463  else
464  CHKERR mField.get_moab().write_file("out_ho_mesh.h5m");
466 }
467 //! [Postprocess results]
468 
469 int main(int argc, char *argv[]) {
470 
471  // Initialisation of MoFEM/PETSc and MOAB data structures
472  const char param_file[] = "param_file.petsc";
473  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
474 
475  auto core_log = logging::core::get();
476  core_log->add_sink(
478  LogManager::setLog("EXAMPLE");
479  MOFEM_LOG_TAG("EXAMPLE", "example");
480 
481  try {
482 
483  //! [Register MoFEM discrete manager in PETSc]
484  DMType dm_name = "DMMOFEM";
485  CHKERR DMRegister_MoFEM(dm_name);
486  //! [Register MoFEM discrete manager in PETSc
487 
488  //! [Create MoAB]
489  moab::Core mb_instance; ///< mesh database
490  moab::Interface &moab = mb_instance; ///< mesh database interface
491  //! [Create MoAB]
492 
493  //! [Create MoFEM]
494  MoFEM::Core core(moab); ///< finite element database
495  MoFEM::Interface &m_field = core; ///< finite element database insterface
496  //! [Create MoFEM]
497 
498  //! [ApproxSphere]
499  ApproxSphere ex(m_field);
500  CHKERR ex.runProblem();
501  //! [ApproxSphere]
502  }
503  CATCH_ERRORS;
504 
506 }
std::string param_file
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()
constexpr double alpha
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
@ 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
#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
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
auto integration_rule
constexpr auto t_kd
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
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:544
boost::ptr_vector< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
boost::ptr_vector< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face 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< 'l', 3 > l
const double T
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
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)
auto createSmartVectorMPI
Create MPI Vector.
CoreTmp< 0 > Core
Definition: Core.hpp:1096
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:41
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]
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::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Basic algebra on fields.
Definition: FieldBlas.hpp:31
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
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.
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:33
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
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.
Postprocess on face.