v0.14.0
Loading...
Searching...
No Matches
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
12using namespace MoFEM;
13
14static char help[] = "...\n\n";
15
17
18constexpr char FIELD_NAME[] = "U";
19constexpr int FIELD_DIM = 1;
20constexpr int SPACE_DIM = 2;
21
22template <int DIM> struct ElementsAndOps {};
23
24template <> struct ElementsAndOps<2> {
26};
27
28template <> struct ElementsAndOps<3> {
30};
31
35
37
38template <int FIELD_DIM> struct ApproxFieldFunction;
39
40template <> 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
50
51struct Example {
52
53 Example(MoFEM::Interface &m_field) : mField(m_field) {}
54
56
57private:
60
62
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
86template <> 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]
147}
148//! [Run programme]
149
150//! [Read mesh]
153
157
159}
160//! [Read mesh]
161
162//! [Set up problem]
165 // Add field
168 constexpr int order = 4;
172}
173//! [Set up problem]
174
175//! [Set integration rule]
178
179 auto rule = [](int, int, int p) -> int { return 2 * p; };
180
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>();
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]
209 auto beta = [](const double, const double, const double) { return 1; };
210 pipeline_mng->getOpDomainLhsPipeline().push_back(
212 pipeline_mng->getOpDomainRhsPipeline().push_back(
215}
216//! [Push operators to pipeline]
217
218//! [Solve]
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]
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
305int main(int argc, char *argv[]) {
306
307 // Initialisation of MoFEM/PETSc and MOAB data structures
308 const char param_file[] = "param_file.petsc";
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 }
334
336}
static Index< 'p', 3 > p
std::string param_file
int main()
Definition: adol-c_atom.cpp:46
static char help[]
constexpr char FIELD_NAME[]
constexpr int SPACE_DIM
constexpr int FIELD_DIM
static const double eps
constexpr char FIELD_NAME[]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
constexpr int FIELD_DIM
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
@ 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_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.
double D
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
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.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
FTensor::Index< 'i', 3 > i
double operator()(const double x, const double y, const double z)
SmartPetscObj< Vec > L2Vec
SmartPetscObj< Vec > resVec
boost::shared_ptr< VectorDouble > approxVals
OpError(boost::shared_ptr< CommonData > &common_data_ptr)
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Example]
Definition: plastic.cpp:226
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode assembleSystem()
MoFEMErrorCode readMesh()
MoFEMErrorCode setIntegrationRules()
static ApproxFieldFunction< FIELD_DIM > approxFunction
boost::shared_ptr< CommonData > commonDataPtr
Definition: integration.cpp:42
MoFEMErrorCode checkResults()
MoFEMErrorCode solveSystem()
MoFEMErrorCode createCommonData()
Example(MoFEM::Interface &m_field)
Simple * simpleInterface
Definition: helmholtz.cpp:41
MoFEMErrorCode runProblem()
MoFEM::Interface & mField
Definition: plastic.cpp:233
MoFEMErrorCode setupProblem()
MoFEMErrorCode outputResults()
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
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.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Get value at integration points for scalar field.
Post post-proc data at points from hash maps.
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 loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition: Simple.cpp:194
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 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.