v0.14.0
approximaton.cpp
Go to the documentation of this file.
1 /**
2  * \file approximation.cpp
3  * \example approximation.cpp
4  *
5  * Using PipelineManager interface calculate the divergence of base functions,
6  * and integral of flux on the boundary. Since the h-div space is used, volume
7  * integral and boundary integral should give the same result.
8  */
9 
10 #include <MoFEM.hpp>
11 
12 using namespace MoFEM;
13 
14 static char help[] = "...\n\n";
15 
16 #include <BasicFiniteElements.hpp>
17 
18 constexpr char FIELD_NAME[] = "U";
19 constexpr int FIELD_DIM = 1;
20 constexpr int SPACE_DIM = 2;
21 
22 template <int DIM> struct ElementsAndOps {};
23 
24 template <> struct ElementsAndOps<2> {
26 };
27 
28 template <> struct ElementsAndOps<3> {
30 };
31 
35 
37 
38 template <int FIELD_DIM> struct ApproxFieldFunction;
39 
40 template <> struct ApproxFieldFunction<1> {
41  double operator()(const double x, const double y, const double z) {
42  return sin(x * 10.) * cos(y * 10.);
43  }
44 };
45 
49  PETSC>::LinearForm<GAUSS>::OpSource<1, FIELD_DIM>;
50 
51 struct Example {
52 
53  Example(MoFEM::Interface &m_field) : mField(m_field) {}
54 
55  MoFEMErrorCode runProblem();
56 
57 private:
58  MoFEM::Interface &mField;
59  Simple *simpleInterface;
60 
62 
63  MoFEMErrorCode readMesh();
64  MoFEMErrorCode setupProblem();
65  MoFEMErrorCode setIntegrationRules();
66  MoFEMErrorCode createCommonData();
67  MoFEMErrorCode boundaryCondition();
68  MoFEMErrorCode assembleSystem();
69  MoFEMErrorCode solveSystem();
70  MoFEMErrorCode outputResults();
71  MoFEMErrorCode checkResults();
72 
73  struct CommonData {
74  boost::shared_ptr<VectorDouble> approxVals;
77  };
78  boost::shared_ptr<CommonData> commonDataPtr;
79 
80  template <int FIELD_DIM> struct OpError;
81 };
82 
85 
86 template <> struct Example::OpError<1> : public DomainEleOp {
87  boost::shared_ptr<CommonData> commonDataPtr;
88  OpError(boost::shared_ptr<CommonData> &common_data_ptr)
89  : DomainEleOp(FIELD_NAME, OPROW), commonDataPtr(common_data_ptr) {}
90  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
92 
93  if (const size_t nb_dofs = data.getIndices().size()) {
94 
95  const int nb_integration_pts = getGaussPts().size2();
96  auto t_w = getFTensor0IntegrationWeight();
97  auto t_val = getFTensor0FromVec(*(commonDataPtr->approxVals));
98  auto t_coords = getFTensor1CoordsAtGaussPts();
99 
100  VectorDouble nf(nb_dofs, false);
101  nf.clear();
102 
104  const double volume = getMeasure();
105 
106  auto t_row_base = data.getFTensor0N();
107  double error = 0;
108  for (int gg = 0; gg != nb_integration_pts; ++gg) {
109 
110  const double alpha = t_w * volume;
111  double diff = t_val - Example::approxFunction(t_coords(0), t_coords(1),
112  t_coords(2));
113  error += alpha * pow(diff, 2);
114 
115  for (size_t r = 0; r != nb_dofs; ++r) {
116  nf[r] += alpha * t_row_base * diff;
117  ++t_row_base;
118  }
119 
120  ++t_w;
121  ++t_val;
122  ++t_coords;
123  }
124 
125  const int index = 0;
126  CHKERR VecSetValue(commonDataPtr->L2Vec, index, error, ADD_VALUES);
127  CHKERR VecSetValues(commonDataPtr->resVec, data, &nf[0], ADD_VALUES);
128  }
129 
131  }
132 };
133 
134 //! [Run programme]
137  CHKERR readMesh();
138  CHKERR setupProblem();
139  CHKERR setIntegrationRules();
140  CHKERR createCommonData();
141  CHKERR boundaryCondition();
142  CHKERR assembleSystem();
143  CHKERR solveSystem();
144  CHKERR outputResults();
145  CHKERR checkResults();
147 }
148 //! [Run programme]
149 
150 //! [Read mesh]
153 
154  CHKERR mField.getInterface(simpleInterface);
155  CHKERR simpleInterface->getOptions();
156  CHKERR simpleInterface->loadFile();
157 
159 }
160 //! [Read mesh]
161 
162 //! [Set up problem]
165  // Add field
166  CHKERR simpleInterface->addDomainField(FIELD_NAME, H1,
168  constexpr int order = 4;
169  CHKERR simpleInterface->setFieldOrder(FIELD_NAME, order);
170  CHKERR simpleInterface->setUp();
172 }
173 //! [Set up problem]
174 
175 //! [Set integration rule]
178 
179  auto rule = [](int, int, int p) -> int { return 2 * p; };
180 
181  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
182  CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule);
183  CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
184 
186 }
187 //! [Set integration rule]
188 
189 //! [Create common data]
192  commonDataPtr = boost::make_shared<CommonData>();
193  commonDataPtr->resVec = createDMVector(simpleInterface->getDM());
194  commonDataPtr->L2Vec = createVectorMPI(
195  mField.get_comm(), (!mField.get_comm_rank()) ? 1 : 0, 1);
196  commonDataPtr->approxVals = boost::make_shared<VectorDouble>();
198 }
199 //! [Create common data]
200 
201 //! [Boundary condition]
203 //! [Boundary condition]
204 
205 //! [Push operators to pipeline]
208  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
209  auto beta = [](const double, const double, const double) { return 1; };
210  pipeline_mng->getOpDomainLhsPipeline().push_back(
211  new OpDomainMass(FIELD_NAME, FIELD_NAME, beta));
212  pipeline_mng->getOpDomainRhsPipeline().push_back(
213  new OpDomainSource(FIELD_NAME, approxFunction));
215 }
216 //! [Push operators to pipeline]
217 
218 //! [Solve]
221  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
222  auto solver = pipeline_mng->createKSP();
223  CHKERR KSPSetFromOptions(solver);
224  CHKERR KSPSetUp(solver);
225 
226  auto dm = simpleInterface->getDM();
227  auto D = createDMVector(dm);
228  auto F = vectorDuplicate(D);
229 
230  CHKERR KSPSolve(solver, F, D);
231  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
232  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
233  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
235 }
236 
237 //! [Solve]
240  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
241  pipeline_mng->getDomainLhsFE().reset();
242  auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
243 
244  auto u_ptr = boost::make_shared<VectorDouble>();
245  post_proc_fe->getOpPtrVector().push_back(
247 
249 
250  post_proc_fe->getOpPtrVector().push_back(
251 
252  new OpPPMap(
253 
254  post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
255 
256  {{FIELD_NAME, u_ptr}},
257 
258  {},
259 
260  {},
261 
262  {})
263 
264  );
265 
266  pipeline_mng->getDomainRhsFE() = post_proc_fe;
267  CHKERR pipeline_mng->loopFiniteElements();
268  CHKERR post_proc_fe->writeFile("out_approx.h5m");
270 }
271 //! [Postprocess results]
272 
273 //! [Check results]
277  pipeline_mng->getDomainLhsFE().reset();
278  pipeline_mng->getDomainRhsFE().reset();
279  pipeline_mng->getOpDomainRhsPipeline().clear();
280  pipeline_mng->getOpDomainRhsPipeline().push_back(
282  pipeline_mng->getOpDomainRhsPipeline().push_back(
284  CHKERR pipeline_mng->loopFiniteElements();
285  CHKERR VecAssemblyBegin(commonDataPtr->L2Vec);
286  CHKERR VecAssemblyEnd(commonDataPtr->L2Vec);
287  CHKERR VecAssemblyBegin(commonDataPtr->resVec);
288  CHKERR VecAssemblyEnd(commonDataPtr->resVec);
289  double nrm2;
290  CHKERR VecNorm(commonDataPtr->resVec, NORM_2, &nrm2);
291  const double *array;
292  CHKERR VecGetArrayRead(commonDataPtr->L2Vec, &array);
293  if (mField.get_comm_rank() == 0)
294  PetscPrintf(PETSC_COMM_SELF, "Error %6.4e Vec norm %6.4e\n", std::sqrt(array[0]),
295  nrm2);
296  CHKERR VecRestoreArrayRead(commonDataPtr->L2Vec, &array);
297  constexpr double eps = 1e-8;
298  if (nrm2 > eps)
299  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
300  "Not converged solution");
302 }
303 //! [Check results]
304 
305 int main(int argc, char *argv[]) {
306 
307  // Initialisation of MoFEM/PETSc and MOAB data structures
308  const char param_file[] = "param_file.petsc";
309  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
310 
311  try {
312 
313  //! [Register MoFEM discrete manager in PETSc]
314  DMType dm_name = "DMMOFEM";
315  CHKERR DMRegister_MoFEM(dm_name);
316  //! [Register MoFEM discrete manager in PETSc
317 
318  //! [Create MoAB]
319  moab::Core mb_instance; ///< mesh database
320  moab::Interface &moab = mb_instance; ///< mesh database interface
321  //! [Create MoAB]
322 
323  //! [Create MoFEM]
324  MoFEM::Core core(moab); ///< finite element database
325  MoFEM::Interface &m_field = core; ///< finite element database insterface
326  //! [Create MoFEM]
327 
328  //! [Example]
329  Example ex(m_field);
330  CHKERR ex.runProblem();
331  //! [Example]
332  }
333  CATCH_ERRORS;
334 
336 }
main
int main(int argc, char *argv[])
[Check results]
Definition: approximaton.cpp:305
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
Example::approxFunction
static ApproxFieldFunction< FIELD_DIM > approxFunction
Definition: approximaton.cpp:61
SPACE_DIM
constexpr int SPACE_DIM
Definition: approximaton.cpp:20
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
Example::checkResults
MoFEMErrorCode checkResults()
[Postprocess results]
Definition: dynamic_first_order_con_law.cpp:1205
MoFEM::PipelineManager::getDomainRhsFE
boost::shared_ptr< FEMethod > & getDomainRhsFE()
Definition: PipelineManager.hpp:405
Example::assembleSystem
MoFEMErrorCode assembleSystem()
[Push operators to pipeline]
Definition: dynamic_first_order_con_law.cpp:647
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
OpDomainSource
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
Definition: approximaton.cpp:49
OpError
Definition: initial_diffusion.cpp:61
Example::Example
Example(MoFEM::Interface &m_field)
Definition: approximaton.cpp:53
Example::CommonData::L2Vec
SmartPetscObj< Vec > L2Vec
Definition: approximaton.cpp:75
FIELD_DIM
constexpr int FIELD_DIM
Definition: approximaton.cpp:19
MoFEM::PipelineManager::loopFiniteElements
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
Definition: PipelineManager.cpp:19
ApproxFieldFunction< 1 >::operator()
double operator()(const double x, const double y, const double z)
Definition: approximaton.cpp:41
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::PipelineManager::getOpDomainRhsPipeline
boost::ptr_deque< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
Definition: PipelineManager.hpp:773
Example::OpError< 1 >::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Definition: approximaton.cpp:90
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
Example::CommonData::approxVals
boost::shared_ptr< VectorDouble > approxVals
Definition: approximaton.cpp:74
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MoFEM.hpp
MoFEM::DMoFEMMeshToLocalVector
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:523
BasicFiniteElements.hpp
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
OpDomainMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
Definition: child_and_parent.cpp:53
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1576
sdf.r
int r
Definition: sdf.py:8
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
OpDomainMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
Definition: approximaton.cpp:47
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition: EntitiesFieldData.hpp:1489
Example::readMesh
MoFEMErrorCode readMesh()
[Run problem]
Definition: dynamic_first_order_con_law.cpp:463
MoFEM::PostProcBrokenMeshInMoab
Definition: PostProcBrokenMeshInMoabBase.hpp:667
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1067
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
Example
[Example]
Definition: plastic.cpp:177
help
static char help[]
Definition: approximaton.cpp:14
double
convert.type
type
Definition: convert.py:64
MoFEM::PipelineManager::FaceEle
MoFEM::FaceElementForcesAndSourcesCore FaceEle
Definition: PipelineManager.hpp:35
Example::OpError< 1 >::OpError
OpError(boost::shared_ptr< CommonData > &common_data_ptr)
Definition: approximaton.cpp:88
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:310
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
Example::solveSystem
MoFEMErrorCode solveSystem()
[Solve]
Definition: dynamic_first_order_con_law.cpp:893
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
MoFEM::PipelineManager::createKSP
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
Definition: PipelineManager.cpp:49
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: EntitiesFieldData.hpp:1201
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
Example::OpError
Definition: approximaton.cpp:80
Example::CommonData
[Example]
Definition: integration.cpp:53
Example::setIntegrationRules
MoFEMErrorCode setIntegrationRules()
[Set up problem]
Definition: plot_base.cpp:203
MoFEM::PipelineManager::setDomainRhsIntegrationRule
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
Definition: PipelineManager.hpp:530
Example::OpError< 1 >::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: approximaton.cpp:87
MoFEM::PipelineManager::getOpDomainLhsPipeline
boost::ptr_deque< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
Definition: PipelineManager.hpp:749
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
BiLinearForm
FTensor::Index< 'i', 3 >
MoFEM::PipelineManager::setDomainLhsIntegrationRule
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Definition: PipelineManager.hpp:503
ElementsAndOps
Definition: child_and_parent.cpp:18
MoFEM::PipelineManager::getDomainLhsFE
boost::shared_ptr< FEMethod > & getDomainLhsFE()
Definition: PipelineManager.hpp:401
DomainEleOp
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
FIELD_NAME
constexpr char FIELD_NAME[]
Definition: approximaton.cpp:18
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:221
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1102
Example::setupProblem
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:225
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
Example::boundaryCondition
MoFEMErrorCode boundaryCondition()
[Set up problem]
Definition: dynamic_first_order_con_law.cpp:536
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_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
EntData
EntitiesFieldData::EntData EntData
Definition: approximaton.cpp:34
MoFEM::createVectorMPI
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
Definition: PetscSmartObj.hpp:202
Example::runProblem
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:213
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
Example::mField
MoFEM::Interface & mField
Definition: plastic.cpp:184
Example::createCommonData
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: plastic.cpp:396
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: child_and_parent.cpp:34
MoFEM::SmartPetscObj< Vec >
OpDomainSource
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
Definition: child_and_parent.cpp:55
Example::outputResults
MoFEMErrorCode outputResults()
[Solve]
Definition: dynamic_first_order_con_law.cpp:1183
Example::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: integration.cpp:40
convert.int
int
Definition: convert.py:64
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
Example::CommonData::resVec
SmartPetscObj< Vec > resVec
Definition: approximaton.cpp:76
ApproxFieldFunction
Definition: child_and_parent.cpp:43
F
@ F
Definition: free_surface.cpp:394
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698