v0.14.0
Loading...
Searching...
No Matches
Classes | Functions | Variables
prism_polynomial_approximation.cpp File Reference
#include <MoFEM.hpp>

Go to the source code of this file.

Classes

struct  ApproxFunction
 
struct  PrismOpCheck
 
struct  PrismOpRhs
 
struct  PrismOpLhs
 
struct  PrismFE
 

Functions

int main (int argc, char *argv[])
 

Variables

static char help [] = "...\n\n"
 
static constexpr int approx_order = 6
 

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)

Definition at line 96 of file prism_polynomial_approximation.cpp.

96 {
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}
@ 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
#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
FTensor::Index< 'n', SPACE_DIM > n
@ 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
double D
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
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
auto createKSP(MPI_Comm comm)
constexpr AssemblyType A
static char help[]
static constexpr int approx_order
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
Deprecated interface functions.
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

Variable Documentation

◆ approx_order

constexpr int approx_order = 6
staticconstexpr

◆ help

char help[] = "...\n\n"
static

Definition at line 12 of file prism_polynomial_approximation.cpp.