v0.9.1
basic_approx.cpp
Go to the documentation of this file.
1 /**
2  * \file basic_approx.cpp
3  * \example basic_approx.cpp
4  *
5  * Using Basic interface calculate the divergence of base functions, and
6  * 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 
35 
36 #include <BaseOps.hpp>
37 
40 
41 struct Example {
42 
43  Example(MoFEM::Interface &m_field) : mField(m_field) {}
44 
45  MoFEMErrorCode runProblem();
46 
47 private:
49 
50  static double approxFunction(const double x, const double y, const double z) {
51  return sin(x * 10.) * cos(y * 10.);
52  }
53  static int integrationRule(int, int, int p_data) { return 2 * p_data; };
54 
55  MoFEMErrorCode setUP();
56  MoFEMErrorCode createCommonData();
57  MoFEMErrorCode bC();
58  MoFEMErrorCode OPs();
59  MoFEMErrorCode kspSolve();
60  MoFEMErrorCode postProcess();
61  MoFEMErrorCode checkResults();
62 
63  struct CommonData {
64  boost::shared_ptr<VectorDouble> approxVals;
67  };
68  boost::shared_ptr<CommonData> commonDataPtr;
69 
70  struct OpError : public DomainEleOp {
71  boost::shared_ptr<CommonData> commonDataPtr;
72  OpError(boost::shared_ptr<CommonData> &common_data_ptr)
73  : DomainEleOp("U", OPROW), commonDataPtr(common_data_ptr) {}
74  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
75  };
76 };
77 
80  CHKERR setUP();
81  CHKERR createCommonData();
82  CHKERR bC();
83  CHKERR OPs();
84  CHKERR kspSolve();
85  CHKERR postProcess();
86  CHKERR checkResults();
88 }
89 
90 //! [Set up problem]
93  Simple *simple = mField.getInterface<Simple>();
94  // Add field
95  CHKERR simple->addDomainField("U", H1, AINSWORTH_BERNSTEIN_BEZIER_BASE, 1);
96  constexpr int order = 4;
97  CHKERR simple->setFieldOrder("U", order);
98  CHKERR simple->setUp();
100 }
101 //! [Set up problem]
102 
103 //! [Create common data]
106  Simple *simple = mField.getInterface<Simple>();
107  commonDataPtr = boost::make_shared<CommonData>();
108  commonDataPtr->resVec = smartCreateDMDVector(simple->getDM());
109  commonDataPtr->L2Vec = createSmartVectorMPI(
110  mField.get_comm(), (!mField.get_comm_rank()) ? 1 : 0, 1);
111  commonDataPtr->approxVals = boost::make_shared<VectorDouble>();
113 }
114 //! [Create common data]
115 
116 //! [Boundary condition]
118 //! [Boundary condition]
119 
120 //! [Push operators to pipeline]
123  Basic *basic = mField.getInterface<Basic>();
124  auto beta = [](const double, const double, const double) { return 1; };
125  basic->getOpDomainLhsPipeline().push_back(new OpDomainMass("U", "U", beta));
126  basic->getOpDomainRhsPipeline().push_back(
128  CHKERR basic->setDomainRhsIntegrationRule(integrationRule);
129  CHKERR basic->setDomainLhsIntegrationRule(integrationRule);
131 }
132 //! [Push operators to pipeline]
133 
134 //! [Solve]
137  Simple *simple = mField.getInterface<Simple>();
138  Basic *basic = mField.getInterface<Basic>();
139  auto solver = basic->createKSP();
140  CHKERR KSPSetFromOptions(solver);
141  CHKERR KSPSetUp(solver);
142 
143  auto dm = simple->getDM();
144  auto D = smartCreateDMDVector(dm);
145  auto F = smartVectorDuplicate(D);
146 
147  CHKERR KSPSolve(solver, F, D);
148  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
149  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
150  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
152 }
153 
154 //! [Solve]
157  Basic *basic = mField.getInterface<Basic>();
158  basic->getDomainLhsFE().reset();
159  auto post_proc_fe = boost::make_shared<PostProcFaceOnRefinedMesh>(mField);
160  post_proc_fe->generateReferenceElementMesh();
161  post_proc_fe->addFieldValuesPostProc("U");
162  basic->getDomainRhsFE() = post_proc_fe;
163  CHKERR basic->loopFiniteElements();
164  CHKERR post_proc_fe->writeFile("out_approx.h5m");
166 }
167 //! [Postprocess results]
168 
169 //! [Solve]
172  Basic *basic = mField.getInterface<Basic>();
173  basic->getDomainLhsFE().reset();
174  basic->getDomainRhsFE().reset();
175  basic->getOpDomainRhsPipeline().clear();
176  basic->getOpDomainRhsPipeline().push_back(
177  new OpCalculateScalarFieldValues("U", commonDataPtr->approxVals));
178  basic->getOpDomainRhsPipeline().push_back(new OpError(commonDataPtr));
179  CHKERR basic->setDomainRhsIntegrationRule(integrationRule);
180  CHKERR basic->loopFiniteElements();
181  CHKERR VecAssemblyBegin(commonDataPtr->L2Vec);
182  CHKERR VecAssemblyEnd(commonDataPtr->L2Vec);
183  CHKERR VecAssemblyBegin(commonDataPtr->resVec);
184  CHKERR VecAssemblyEnd(commonDataPtr->resVec);
185  double nrm2;
186  CHKERR VecNorm(commonDataPtr->resVec, NORM_2, &nrm2);
187  const double *array;
188  CHKERR VecGetArrayRead(commonDataPtr->L2Vec, &array);
189  if (mField.get_comm_rank() == 0)
190  PetscPrintf(PETSC_COMM_SELF, "Error %6.4e Vec norm %6.4e\n", sqrt(array[0]),
191  nrm2);
192  CHKERR VecRestoreArrayRead(commonDataPtr->L2Vec, &array);
193  constexpr double eps = 1e-8;
194  if (nrm2 > eps)
195  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
196  "Not converged solution");
198 }
199 //! [Solver]
200 
201 int main(int argc, char *argv[]) {
202 
203  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
204 
205  try {
206 
207  //! [Register MoFEM discrete manager in PETSc]
208  DMType dm_name = "DMMOFEM";
209  CHKERR DMRegister_MoFEM(dm_name);
210  //! [Register MoFEM discrete manager in PETSc
211 
212  //! [Create MoAB]
213  moab::Core mb_instance; ///< mesh database
214  moab::Interface &moab = mb_instance; ///< mesh database interface
215  //! [Create MoAB]
216 
217  //! [Create MoFEM]
218  MoFEM::Core core(moab); ///< finite element database
219  MoFEM::Interface &m_field = core; ///< finite element database insterface
220  //! [Create MoFEM]
221 
222  //! [Load mesh]
223  Simple *simple = m_field.getInterface<Simple>();
224  CHKERR simple->getOptions();
225  CHKERR simple->loadFile("");
226  //! [Load mesh]
227 
228  //! [Example]
229  Example ex(m_field);
230  CHKERR ex.runProblem();
231  //! [Example]
232  }
233  CATCH_ERRORS;
234 
236 }
237 
239  EntData &data) {
241 
242  if (const size_t nb_dofs = data.getIndices().size()) {
243 
244  const int nb_integration_pts = getGaussPts().size2();
245  auto t_w = getFTensor0IntegrationWeight();
246  auto t_val = getFTensor0FromVec(*(commonDataPtr->approxVals));
247  auto t_coords = getFTensor1CoordsAtGaussPts();
248 
249  VectorDouble nf(nb_dofs, false);
250  nf.clear();
251 
253  const double volume = getMeasure();
254 
255  auto t_row_base = data.getFTensor0N();
256  double error = 0;
257  for (int gg = 0; gg != nb_integration_pts; ++gg) {
258 
259  const double alpha = t_w * volume;
260  double diff = t_val - Example::approxFunction(t_coords(0), t_coords(1),
261  t_coords(2));
262  error += alpha * pow(diff, 2);
263 
264  for (size_t r = 0; r != nb_dofs; ++r) {
265  nf[r] += alpha * t_row_base * diff;
266  ++t_row_base;
267  }
268 
269  ++t_w;
270  ++t_val;
271  ++t_coords;
272  }
273 
274  const int index = 0;
275  CHKERR VecSetValue(commonDataPtr->L2Vec, index, error, ADD_VALUES);
276  CHKERR VecSetValues(commonDataPtr->resVec, data, &nf[0], ADD_VALUES);
277  }
278 
280 }
boost::shared_ptr< CommonData > commonDataPtr
SmartPetscObj< Vec > resVec
static FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:141
static char help[]
Deprecated interface functions.
MoFEMErrorCode createCommonData()
[Set up problem]
MoFEMErrorCode kspSolve()
[Push operators to pipeline]
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
static double approxFunction(const double x, const double y, const double z)
[Grad residual operator]
Definition: BaseOps.hpp:239
MoFEM::Interface & mField
boost::ptr_vector< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
Definition: Basic.hpp:466
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:431
DataForcesAndSourcesCore::EntData EntData
boost::shared_ptr< FEMethod > & getDomainLhsFE()
Definition: Basic.hpp:272
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
Example(MoFEM::Interface &m_field)
boost::ptr_vector< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
Definition: Basic.hpp:442
Core (interface) class.
Definition: Core.hpp:50
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
Definition: AuxPETSc.hpp:238
OpError(boost::shared_ptr< CommonData > &common_data_ptr)
SmartPetscObj< Vec > L2Vec
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
Simple interface for fast problem set-up.
Definition: Simple.hpp:36
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
Definition: Basic.cpp:68
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:51
auto createSmartVectorMPI
Create MPI Vector.
Definition: AuxPETSc.hpp:224
const VectorInt & getIndices() const
Get global indices of dofs on entity.
MoFEMErrorCode runProblem()
[Run problem]
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
MoFEMErrorCode OPs()
[Boundary condition]
Basic interface.
Definition: Basic.hpp:36
OpTools< DomainEleOp >::OpSource< 2 > OpDomainSource
MoFEMErrorCode checkResults()
[Postprocess results]
MoFEMErrorCode bC()
[Create common data]
constexpr int order
[Source operator]
Definition: BaseOps.hpp:97
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:48
static moab::Error error
static int integrationRule(int, int, int p_data)
MoFEMErrorCode VecSetValues(Vec V, const DataForcesAndSourcesCore::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
Definition: Basic.hpp:274
OpTools< DomainEleOp >::OpMass OpDomainMass
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
Definition: Basic.cpp:36
#define CHKERR
Inline error check.
Definition: definitions.h:601
MoFEMErrorCode setUP()
[Set up problem]
boost::shared_ptr< CommonData > commonDataPtr
boost::shared_ptr< VectorDouble > approxVals
MoFEMErrorCode postProcess()
[Solve]
int main(int argc, char *argv[])
[Solver]
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1791
Data on single entity (This is passed as argument to DataOperator::doWork)
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Definition: Basic.hpp:285
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
Definition: Basic.hpp:311
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:438
continuous field
Definition: definitions.h:176
static const double eps
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72
auto smartCreateDMDVector
Get smart vector from DM.
Definition: DMMoFEM.hpp:907
Get value at integration points for scalar field.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:61
Face finite elementUser is implementing own operator at Gauss point level, by own object derived from...
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:26