v0.14.0
Loading...
Searching...
No Matches
prism_polynomial_approximation.cpp
Go to the documentation of this file.
1/** \file prism_polynomial_approximation.cpp
2 \example prism_polynomial_approximation.cpp
3 \brief Checking approximation functions on prism
4*/
5
6
7
8#include <MoFEM.hpp>
9
10using namespace MoFEM;
11
12static char help[] = "...\n\n";
13
14static constexpr int approx_order = 6;
15
17 static inline double fun(double x, double y, double z) {
18 double r = 1;
19 for (int o = 1; o <= approx_order; ++o) {
20 for (int i = 0; i <= o; ++i) {
21 for (int j = 0; j <= (o - i); ++j) {
22 int k = o - i - j;
23 if (k >= 0) {
24 r += pow(x, i) * pow(y, j) * pow(z, k);
25 }
26 }
27 }
28 }
29 return r;
30 }
31
32 static inline VectorDouble3 diff_fun(double x, double y, double z) {
33 VectorDouble3 r(3);
34 r.clear();
35 for (int o = 1; o <= approx_order; ++o) {
36 for (int i = 0; i <= o; ++i) {
37 for (int j = 0; j <= (o - i); ++j) {
38 int k = o - i - j;
39 if (k >= 0) {
40 r[0] += i > 0 ? i * pow(x, i - 1) * pow(y, j) * pow(z, k) : 0;
41 r[1] += j > 0 ? j * pow(x, i) * pow(y, j - 1) * pow(z, k) : 0;
42 r[2] += k > 0 ? k * pow(x, i) * pow(y, j) * pow(z, k - 1) : 0;
43 }
44 }
45 }
46 }
47 return r;
48 }
49};
50
53
54 PrismOpCheck(boost::shared_ptr<VectorDouble> &field_vals,
55 boost::shared_ptr<MatrixDouble> &diff_field_vals);
56 MoFEMErrorCode doWork(int side, EntityType type,
58
59private:
60 boost::shared_ptr<VectorDouble> fieldVals;
61 boost::shared_ptr<MatrixDouble> diffFieldVals;
62};
63
66
68 MoFEMErrorCode doWork(int side, EntityType type,
70
71private:
73};
74
77
79 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
80 EntityType col_type,
83
84private:
86};
87
89
91 FatPrismElementForcesAndSourcesCore;
94};
95
96int main(int argc, char *argv[]) {
97
98 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
99
100 try {
101
102 moab::Core mb_instance;
103 moab::Interface &moab = mb_instance;
104
105 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
106 auto moab_comm_wrap =
107 boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD, false);
108 if (pcomm == NULL)
109 pcomm = new ParallelComm(&moab, moab_comm_wrap->get_comm());
110
111 std::array<double, 18> one_prism_coords = {0, 0, 0, 1, 0, 0, 0, 1, 0,
112 0, 0, 1, 1, 0, 1, 0, 1, 1};
113 std::array<EntityHandle, 6> one_prism_nodes;
114 for (int n = 0; n != 6; ++n)
115 CHKERR moab.create_vertex(&one_prism_coords[3 * n], one_prism_nodes[n]);
116 EntityHandle one_prism;
117 CHKERR moab.create_element(MBPRISM, one_prism_nodes.data(), 6, one_prism);
118 Range one_prism_range;
119 one_prism_range.insert(one_prism);
120 Range one_prism_adj_ents;
121 for (int d = 1; d != 3; ++d)
122 CHKERR moab.get_adjacencies(one_prism_range, d, true, one_prism_adj_ents,
123 moab::Interface::UNION);
124
125 MoFEM::Core core(moab);
126 MoFEM::Interface &m_field = core;
127
128 PrismsFromSurfaceInterface *prisms_from_surface_interface;
129 CHKERR m_field.getInterface(prisms_from_surface_interface);
130 BitRefLevel bit_level0;
131 bit_level0.set(0);
132 CHKERR m_field.getInterface<BitRefManager>()->setEntitiesBitRefLevel(
133 one_prism_range, bit_level0);
134 CHKERR prisms_from_surface_interface->seedPrismsEntities(one_prism_range,
135 bit_level0);
136
137 // Fields
138 CHKERR m_field.add_field("FIELD1", H1, AINSWORTH_LEGENDRE_BASE, 1);
139 CHKERR m_field.add_ents_to_field_by_type(0, MBPRISM, "FIELD1");
140
141 CHKERR m_field.set_field_order(0, MBVERTEX, "FIELD1", 1);
142 CHKERR m_field.set_field_order(0, MBEDGE, "FIELD1", approx_order);
143 CHKERR m_field.set_field_order(0, MBTRI, "FIELD1", approx_order);
144 CHKERR m_field.set_field_order(0, MBQUAD, "FIELD1", approx_order);
145 CHKERR m_field.set_field_order(0, MBPRISM, "FIELD1", approx_order);
146 CHKERR m_field.build_fields();
147
148 // FE
149 CHKERR m_field.add_finite_element("PRISM");
150
151 // Define rows/cols and element data
152 CHKERR m_field.modify_finite_element_add_field_row("PRISM", "FIELD1");
153 CHKERR m_field.modify_finite_element_add_field_col("PRISM", "FIELD1");
154 CHKERR m_field.modify_finite_element_add_field_data("PRISM", "FIELD1");
155 CHKERR m_field.add_ents_to_finite_element_by_type(0, MBPRISM, "PRISM");
156
157 // build finite elemnts
159 // //build adjacencies
160 CHKERR m_field.build_adjacencies(bit_level0);
161
162 // Problem
163 CHKERR m_field.add_problem("TEST_PROBLEM");
164
165 // set finite elements for problem
166 CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM", "PRISM");
167 // set refinement level for problem
168 CHKERR m_field.modify_problem_ref_level_add_bit("TEST_PROBLEM", bit_level0);
169
170 // build problem
171 ProblemsManager *prb_mng_ptr;
172 CHKERR m_field.getInterface(prb_mng_ptr);
173 CHKERR prb_mng_ptr->buildProblem("TEST_PROBLEM", true);
174 // partition
175 CHKERR prb_mng_ptr->partitionSimpleProblem("TEST_PROBLEM");
176 CHKERR prb_mng_ptr->partitionFiniteElements("TEST_PROBLEM");
177 // what are ghost nodes, see Petsc Manual
178 CHKERR prb_mng_ptr->partitionGhostDofs("TEST_PROBLEM");
179
180 // Create matrices
183 ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("TEST_PROBLEM", A);
185 CHKERR m_field.getInterface<VecManager>()->vecCreateGhost("TEST_PROBLEM",
186 ROW, F);
188 CHKERR m_field.getInterface<VecManager>()->vecCreateGhost("TEST_PROBLEM",
189 COL, D);
190
191 auto assemble_matrices_and_vectors = [&]() {
193 PrismFE fe(m_field);
194 MatrixDouble inv_jac;
195 fe.getOpPtrVector().push_back(
197 fe.getOpPtrVector().push_back(
199 fe.getOpPtrVector().push_back(
201 fe.getOpPtrVector().push_back(new PrismOpRhs(F));
202 fe.getOpPtrVector().push_back(new PrismOpLhs(A));
203 CHKERR VecZeroEntries(F);
204 CHKERR MatZeroEntries(A);
205 CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "PRISM", fe);
206 CHKERR VecAssemblyBegin(F);
207 CHKERR VecAssemblyEnd(F);
208 CHKERR MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
209 CHKERR MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
211 };
212
213 auto solve_problem = [&] {
215 auto solver = createKSP(PETSC_COMM_WORLD);
216 CHKERR KSPSetOperators(solver, A, A);
217 CHKERR KSPSetFromOptions(solver);
218 CHKERR KSPSetUp(solver);
219 CHKERR KSPSolve(solver, F, D);
220 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
221 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
222 CHKERR m_field.getInterface<VecManager>()->setLocalGhostVector(
223 "TEST_PROBLEM", COL, D, INSERT_VALUES, SCATTER_REVERSE);
225 };
226
227 auto check_solution = [&] {
229 PrismFE fe(m_field);
230 boost::shared_ptr<VectorDouble> field_vals_ptr(new VectorDouble());
231 boost::shared_ptr<MatrixDouble> diff_field_vals_ptr(new MatrixDouble());
232 MatrixDouble inv_jac;
233 fe.getOpPtrVector().push_back(
234 new OpCalculateInvJacForFatPrism(inv_jac));
235 fe.getOpPtrVector().push_back(
236 new OpSetInvJacH1ForFatPrism(inv_jac));
237 fe.getOpPtrVector().push_back(
238 new OpCalculateScalarFieldValues("FIELD1", field_vals_ptr));
239 fe.getOpPtrVector().push_back(
240 new OpCalculateScalarFieldGradient<3>("FIELD1", diff_field_vals_ptr));
241 fe.getOpPtrVector().push_back(
242 new PrismOpCheck(field_vals_ptr, diff_field_vals_ptr));
243 CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "PRISM", fe);
245 };
246
247 CHKERR assemble_matrices_and_vectors();
248 CHKERR solve_problem();
249 CHKERR check_solution();
250 }
252
254 return 0;
255}
256
257PrismOpCheck::PrismOpCheck(boost::shared_ptr<VectorDouble> &field_vals,
258 boost::shared_ptr<MatrixDouble> &diff_field_vals)
260 "FIELD1", "FIELD1", ForcesAndSourcesCore::UserDataOperator::OPROW),
261 fieldVals(field_vals), diffFieldVals(diff_field_vals) {}
262
266 if (type == MBVERTEX) {
267 const int nb_gauss_pts = data.getN().size2();
268 auto t_coords = getFTensor1CoordsAtGaussPts();
269 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
270 double f = ApproxFunction::fun(t_coords(0), t_coords(1), t_coords(2));
271 VectorDouble3 diff_f =
272 ApproxFunction::diff_fun(t_coords(0), t_coords(1), t_coords(2));
273
274 std::cout << f - (*fieldVals)[gg] << " : ";
275 for (auto d : {0, 1, 2})
276 std::cout << diff_f[d] - (*diffFieldVals)(d, gg) << " ";
277 std::cout << std::endl;
278
279 constexpr double eps = 1e-6;
280 if (std::abs(f - (*fieldVals)[gg]) > eps ||
281 !std::isnormal((*fieldVals)[gg]))
282 SETERRQ3(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
283 "Wrong value %6.4e != %6.4e (%6.4e)", f, (*fieldVals)[gg],
284 f - (*fieldVals)[gg]);
285
286 for (auto d : {0, 1, 2})
287 if (std::abs(diff_f[d] - (*diffFieldVals)(d, gg)) > eps ||
288 !std::isnormal((*diffFieldVals)(d, gg)))
289 SETERRQ3(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
290 "Wrong diff value %6.4e != %6.4e (%6.4e)", diff_f[d],
291 (*diffFieldVals)(d, gg),
292 diff_f[d] - (*diffFieldVals)(d, gg));
293
294 ++t_coords;
295 }
296 }
298}
299
302 "FIELD1", "FIELD1", ForcesAndSourcesCore::UserDataOperator::OPROW),
303 F(f) {}
304
308 const int nb_dofs = data.getN().size2();
309 if (nb_dofs) {
310 const int nb_gauss_pts = data.getN().size1();
311 VectorDouble nf(nb_dofs);
312 nf.clear();
313 auto t_base = data.getFTensor0N();
314 auto t_coords = getFTensor1CoordsAtGaussPts();
315 auto t_w = getFTensor0IntegrationWeight();
316 double vol = getVolume();
317 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
318 double f = ApproxFunction::fun(t_coords(0), t_coords(1), t_coords(2));
319 double v = t_w * vol * f;
320 double *val = &*nf.begin();
321 for (int bb = 0; bb != nb_dofs; ++bb) {
322 *val += v * t_base;
323 ++t_base;
324 ++val;
325 }
326 ++t_coords;
327 ++t_w;
328 }
329 CHKERR VecSetValues(F, data, &*nf.data().begin(), ADD_VALUES);
330 }
332}
333
336 "FIELD1", "FIELD1", ForcesAndSourcesCore::UserDataOperator::OPROWCOL),
337 A(a) {
338 // FIXME: Can be symmetric, is not for simplicity
339 sYmm = false;
340}
341
342MoFEMErrorCode PrismOpLhs::doWork(int row_side, int col_side,
343 EntityType row_type, EntityType col_type,
345 EntitiesFieldData::EntData &col_data) {
347 const int row_nb_dofs = row_data.getN().size2();
348 const int col_nb_dofs = col_data.getN().size2();
349 if (row_nb_dofs && col_nb_dofs) {
350 const int nb_gauss_pts = row_data.getN().size1();
351 MatrixDouble m(row_nb_dofs, col_nb_dofs);
352 m.clear();
353 auto t_w = getFTensor0IntegrationWeight();
354 double vol = getVolume();
355 double *row_base_ptr = &*row_data.getN().data().begin();
356 double *col_base_ptr = &*col_data.getN().data().begin();
357 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
358
359 double v = t_w * vol;
360 cblas_dger(CblasRowMajor, row_nb_dofs, col_nb_dofs, v, row_base_ptr, 1,
361 col_base_ptr, 1, &*m.data().begin(), col_nb_dofs);
362
363 row_base_ptr += row_nb_dofs;
364 col_base_ptr += col_nb_dofs;
365 ++t_w;
366 }
367 CHKERR MatSetValues(A, row_data, col_data, &*m.data().begin(), ADD_VALUES);
368 }
370}
371
372int PrismFE::getRuleTrianglesOnly(int order) { return 2 * (order + 1); };
373int PrismFE::getRuleThroughThickness(int order) { return 2 * (order + 1); };
static Index< 'o', 3 > o
ForcesAndSourcesCore::UserDataOperator UserDataOperator
int main()
Definition: adol-c_atom.cpp:46
constexpr double a
static const double eps
@ COL
Definition: definitions.h:123
@ ROW
Definition: definitions.h:123
#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 MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
#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
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
@ F
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
virtual MoFEMErrorCode loop_finite_elements(const std::string problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
MoFEMErrorCode partitionGhostDofs(const std::string name, int verb=VERBOSE)
determine ghost nodes
MoFEMErrorCode partitionSimpleProblem(const std::string name, int verb=VERBOSE)
partition problem dofs
MoFEMErrorCode buildProblem(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures
MoFEMErrorCode partitionFiniteElements(const std::string name, bool part_from_moab=false, int low_proc=-1, int hi_proc=-1, int verb=VERBOSE)
partition finite elements
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
virtual MoFEMErrorCode modify_problem_add_finite_element(const std::string name_problem, const std::string &fe_name)=0
add finite element to problem, this add entities assigned to finite element to a particular problem
FTensor::Index< 'i', SPACE_DIM > i
double D
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:92
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
auto createKSP(MPI_Comm comm)
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
constexpr AssemblyType A
static char help[]
static constexpr int approx_order
static constexpr int approx_order
static VectorDouble3 diff_fun(double x, double y, double z)
static double fun(double x, double y, double z)
Managing BitRefLevels.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode add_field(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_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
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
bool sYmm
If true assume that matrix is symmetric structure.
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.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
auto getFTensor0IntegrationWeight()
Get integration weights.
structure to get information form mofem into EntitiesFieldData
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Matrix manager is used to build and partition problems.
Calculate inverse of jacobian for face element.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Operator for fat prism element updating integration weights in the volume.
Transform local reference derivatives of shape functions to global derivatives.
MoFEMErrorCode seedPrismsEntities(Range &prisms, const BitRefLevel &bit, int verb=-1)
Seed prism entities by bit level.
Problem manager is used to build and partition problems.
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:23
int getRuleThroughThickness(int order)
int getRuleTrianglesOnly(int order)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< MatrixDouble > diffFieldVals
boost::shared_ptr< VectorDouble > fieldVals
PrismOpCheck(boost::shared_ptr< VectorDouble > &field_vals, boost::shared_ptr< MatrixDouble > &diff_field_vals)
SmartPetscObj< Mat > A
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
PrismOpLhs(SmartPetscObj< Mat > &a)
SmartPetscObj< Vec > F
PrismOpRhs(SmartPetscObj< Vec > &f)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.