v0.13.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 /* This file is part of MoFEM.
11  * MoFEM is free software: you can redistribute it and/or modify it under
12  * the terms of the GNU Lesser General Public License as published by the
13  * Free Software Foundation, either version 3 of the License, or (at your
14  * option) any later version.
15  *
16  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
17  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
19  * License for more details.
20  *
21  * You should have received a copy of the GNU Lesser General Public
22  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
23 
24 #include <MoFEM.hpp>
25 
26 using namespace MoFEM;
27 
28 static char help[] = "...\n\n";
29 
30 #include <BasicFiniteElements.hpp>
31 
32 constexpr char FIELD_NAME[] = "U";
33 constexpr int FIELD_DIM = 1;
34 constexpr int SPACE_DIM = 2;
35 
36 template <int DIM> struct ElementsAndOps {};
37 
38 template <> struct ElementsAndOps<2> {
42 };
43 
44 template <> struct ElementsAndOps<3> {
48 };
49 
53 
54 template <int FIELD_DIM> struct ApproxFieldFunction;
55 
56 template <> struct ApproxFieldFunction<1> {
57  double operator()(const double x, const double y, const double z) {
58  return sin(x * 10.) * cos(y * 10.);
59  }
60 };
61 
66 
67 struct Example {
68 
69  Example(MoFEM::Interface &m_field) : mField(m_field) {}
70 
72 
73 private:
74  MoFEM::Interface &mField;
75  Simple *simpleInterface;
76 
78 
88 
89  struct CommonData {
90  boost::shared_ptr<VectorDouble> approxVals;
93  };
94  boost::shared_ptr<CommonData> commonDataPtr;
95 
96  template <int FIELD_DIM> struct OpError;
97 };
98 
101 
102 template <> struct Example::OpError<1> : public DomainEleOp {
103  boost::shared_ptr<CommonData> commonDataPtr;
104  OpError(boost::shared_ptr<CommonData> &common_data_ptr)
105  : DomainEleOp(FIELD_NAME, OPROW), commonDataPtr(common_data_ptr) {}
106  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
108 
109  if (const size_t nb_dofs = data.getIndices().size()) {
110 
111  const int nb_integration_pts = getGaussPts().size2();
112  auto t_w = getFTensor0IntegrationWeight();
113  auto t_val = getFTensor0FromVec(*(commonDataPtr->approxVals));
114  auto t_coords = getFTensor1CoordsAtGaussPts();
115 
116  VectorDouble nf(nb_dofs, false);
117  nf.clear();
118 
120  const double volume = getMeasure();
121 
122  auto t_row_base = data.getFTensor0N();
123  double error = 0;
124  for (int gg = 0; gg != nb_integration_pts; ++gg) {
125 
126  const double alpha = t_w * volume;
127  double diff = t_val - Example::approxFunction(t_coords(0), t_coords(1),
128  t_coords(2));
129  error += alpha * pow(diff, 2);
130 
131  for (size_t r = 0; r != nb_dofs; ++r) {
132  nf[r] += alpha * t_row_base * diff;
133  ++t_row_base;
134  }
135 
136  ++t_w;
137  ++t_val;
138  ++t_coords;
139  }
140 
141  const int index = 0;
142  CHKERR VecSetValue(commonDataPtr->L2Vec, index, error, ADD_VALUES);
143  CHKERR VecSetValues(commonDataPtr->resVec, data, &nf[0], ADD_VALUES);
144  }
145 
147  }
148 };
149 
150 //! [Run programme]
153  CHKERR readMesh();
154  CHKERR setupProblem();
155  CHKERR setIntegrationRules();
156  CHKERR createCommonData();
157  CHKERR boundaryCondition();
158  CHKERR assembleSystem();
159  CHKERR solveSystem();
160  CHKERR outputResults();
161  CHKERR checkResults();
163 }
164 //! [Run programme]
165 
166 //! [Read mesh]
169 
170  CHKERR mField.getInterface(simpleInterface);
171  CHKERR simpleInterface->getOptions();
172  CHKERR simpleInterface->loadFile();
173 
175 }
176 //! [Read mesh]
177 
178 //! [Set up problem]
181  // Add field
182  CHKERR simpleInterface->addDomainField(FIELD_NAME, H1,
184  constexpr int order = 4;
185  CHKERR simpleInterface->setFieldOrder(FIELD_NAME, order);
186  CHKERR simpleInterface->setUp();
188 }
189 //! [Set up problem]
190 
191 //! [Set integration rule]
194 
195  auto rule = [](int, int, int p) -> int { return 2 * p; };
196 
197  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
198  CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule);
199  CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
200 
202 }
203 //! [Set integration rule]
204 
205 //! [Create common data]
208  commonDataPtr = boost::make_shared<CommonData>();
209  commonDataPtr->resVec = smartCreateDMVector(simpleInterface->getDM());
210  commonDataPtr->L2Vec = createSmartVectorMPI(
211  mField.get_comm(), (!mField.get_comm_rank()) ? 1 : 0, 1);
212  commonDataPtr->approxVals = boost::make_shared<VectorDouble>();
214 }
215 //! [Create common data]
216 
217 //! [Boundary condition]
219 //! [Boundary condition]
220 
221 //! [Push operators to pipeline]
224  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
225  auto beta = [](const double, const double, const double) { return 1; };
226  pipeline_mng->getOpDomainLhsPipeline().push_back(
227  new OpDomainMass(FIELD_NAME, FIELD_NAME, beta));
228  pipeline_mng->getOpDomainRhsPipeline().push_back(
229  new OpDomainSource(FIELD_NAME, approxFunction));
231 }
232 //! [Push operators to pipeline]
233 
234 //! [Solve]
237  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
238  auto solver = pipeline_mng->createKSP();
239  CHKERR KSPSetFromOptions(solver);
240  CHKERR KSPSetUp(solver);
241 
242  auto dm = simpleInterface->getDM();
243  auto D = smartCreateDMVector(dm);
244  auto F = smartVectorDuplicate(D);
245 
246  CHKERR KSPSolve(solver, F, D);
247  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
248  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
249  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
251 }
252 
253 //! [Solve]
256  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
257  pipeline_mng->getDomainLhsFE().reset();
258  auto post_proc_fe = boost::make_shared<PostProcFaceOnRefinedMesh>(mField);
259  post_proc_fe->generateReferenceElementMesh();
260  post_proc_fe->addFieldValuesPostProc(FIELD_NAME);
261  pipeline_mng->getDomainRhsFE() = post_proc_fe;
262  CHKERR pipeline_mng->loopFiniteElements();
263  CHKERR post_proc_fe->writeFile("out_approx.h5m");
265 }
266 //! [Postprocess results]
267 
268 //! [Check results]
271  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
272  pipeline_mng->getDomainLhsFE().reset();
273  pipeline_mng->getDomainRhsFE().reset();
274  pipeline_mng->getOpDomainRhsPipeline().clear();
275  pipeline_mng->getOpDomainRhsPipeline().push_back(
276  new OpCalculateScalarFieldValues(FIELD_NAME, commonDataPtr->approxVals));
277  pipeline_mng->getOpDomainRhsPipeline().push_back(
278  new OpError<FIELD_DIM>(commonDataPtr));
279  CHKERR pipeline_mng->loopFiniteElements();
280  CHKERR VecAssemblyBegin(commonDataPtr->L2Vec);
281  CHKERR VecAssemblyEnd(commonDataPtr->L2Vec);
282  CHKERR VecAssemblyBegin(commonDataPtr->resVec);
283  CHKERR VecAssemblyEnd(commonDataPtr->resVec);
284  double nrm2;
285  CHKERR VecNorm(commonDataPtr->resVec, NORM_2, &nrm2);
286  const double *array;
287  CHKERR VecGetArrayRead(commonDataPtr->L2Vec, &array);
288  if (mField.get_comm_rank() == 0)
289  PetscPrintf(PETSC_COMM_SELF, "Error %6.4e Vec norm %6.4e\n", std::sqrt(array[0]),
290  nrm2);
291  CHKERR VecRestoreArrayRead(commonDataPtr->L2Vec, &array);
292  constexpr double eps = 1e-8;
293  if (nrm2 > eps)
294  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
295  "Not converged solution");
297 }
298 //! [Check results]
299 
300 int main(int argc, char *argv[]) {
301 
302  // Initialisation of MoFEM/PETSc and MOAB data structures
303  const char param_file[] = "param_file.petsc";
304  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
305 
306  try {
307 
308  //! [Register MoFEM discrete manager in PETSc]
309  DMType dm_name = "DMMOFEM";
310  CHKERR DMRegister_MoFEM(dm_name);
311  //! [Register MoFEM discrete manager in PETSc
312 
313  //! [Create MoAB]
314  moab::Core mb_instance; ///< mesh database
315  moab::Interface &moab = mb_instance; ///< mesh database interface
316  //! [Create MoAB]
317 
318  //! [Create MoFEM]
319  MoFEM::Core core(moab); ///< finite element database
320  MoFEM::Interface &m_field = core; ///< finite element database insterface
321  //! [Create MoFEM]
322 
323  //! [Example]
324  Example ex(m_field);
325  CHKERR ex.runProblem();
326  //! [Example]
327  }
328  CATCH_ERRORS;
329 
331 }
static Index< 'p', 3 > p
std::string param_file
int main(int argc, char *argv[])
[Check results]
EntitiesFieldData::EntData EntData
static char help[]
constexpr char FIELD_NAME[]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
constexpr int SPACE_DIM
constexpr int FIELD_DIM
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
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
constexpr double alpha
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ 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
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
#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
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
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.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
FTensor::Index< 'i', SPACE_DIM > i
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
auto createSmartVectorMPI
Create MPI Vector.
CoreTmp< 0 > Core
Definition: Core.hpp:1096
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:149
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
const double D
diffusivity
const double r
rate factor
double operator()(const double x, const double y, const double z)
PipelineManager::FaceEle DomainEle
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)
Operator for linear form, usually to calculate values on right hand side.
[Example]
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode assembleSystem()
MoFEMErrorCode readMesh()
MoFEMErrorCode setIntegrationRules()
static ApproxFieldFunction< FIELD_DIM > approxFunction
MoFEMErrorCode checkResults()
MoFEMErrorCode solveSystem()
MoFEMErrorCode createCommonData()
Example(MoFEM::Interface &m_field)
MoFEMErrorCode runProblem()
MoFEMErrorCode setupProblem()
MoFEMErrorCode outputResults()
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.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Get value at integration points for scalar field.
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:33
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
Postprocess on face.
Post processing.