v0.13.1
test_cache_on_entities.cpp
Go to the documentation of this file.
1/**
2 * \file test_cache_on_entities.cpp
3 * \example test_cache_on_entities.cpp
4 *
5 * Testing entities cache. Entities cache is data abstraction enabling user to
6 * pass data between operators, finite elements, and problems, in convenient
7 * way.
8 *
9 * It can be used to store indices, or history variables, and any other data,
10 * which I can not imagine at that moment of time, and how to make use of it
11 * depends on user creativity and imagination.
12 *
13 */
14
15
16
17#include <MoFEM.hpp>
18using namespace MoFEM;
19
20static char help[] = "...\n\n";
21
26
30};
31
32/**
33 * @brief Operator set cache stored data, in this example, global indices, but
34 * it can be any structure
35 *
36 */
37struct OpVolumeSet : public VolOp {
38 OpVolumeSet(const std::string &field_name) : VolOp(field_name, OPROW) {}
41
43 MOFEM_LOG_CHANNEL("SYNC");
44 MOFEM_LOG_TAG("SYNC", "OpVolumeSet");
45
46 // Clear data when start process element
47 if (type == MBVERTEX)
48 entsIndices.clear();
49
50 // Get field entity pointer
51 auto field_ents = data.getFieldEntities();
52 if (auto e_ptr = field_ents[0]) {
53 // Add indices to global storage
54 entsIndices.push_back(boost::make_shared<MyStorage>(data.getIndices()));
55 // Store pointer to data on entity
56 e_ptr->getWeakStoragePtr() = entsIndices.back();
57
58 MOFEM_LOG("SYNC", Sev::inform) << "Set " << e_ptr->getEntTypeName() << " "
59 << side << " : " << entsIndices.size();
60
61 // Check if all works
62 if (auto ptr = e_ptr->getSharedStoragePtr<EntityStorage>()) {
63
64 if (auto cast_ptr = boost::dynamic_pointer_cast<MyStorage>(ptr))
65 MOFEM_LOG("SYNC", Sev::noisy) << "Cast works";
66 else
67 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Cast not works");
68
69 } else {
70 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
71 "Pointer to stored data is expected to be set");
72 }
73
74 if (auto ptr = e_ptr->getSharedStoragePtr<MyStorage>()) {
75 MOFEM_LOG("SYNC", Sev::verbose)
76 << data.getIndices() << " : " << ptr->globalIndices;
77 } else {
78 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
79 "Pointer to stored data is expected to be set");
80 }
81
82 } else {
83 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
84 "Pointer to entity should be set by finite element");
85 }
86
88 }
89
90 using SharedVecInt = boost::shared_ptr<MyStorage>;
91 static std::vector<SharedVecInt>
92 entsIndices; ///< This is global static storage
93};
94
95std::vector<OpVolumeSet::SharedVecInt> OpVolumeSet::entsIndices;
96
97/**
98 * @brief Test if cached data can be accessed, and check consistency of data
99 *
100 */
101struct OpVolumeTest : public VolOp {
102 OpVolumeTest(const std::string &field_name) : VolOp(field_name, OPROW) {}
105
107 MOFEM_LOG_CHANNEL("SYNC");
108 MOFEM_LOG_TAG("SYNC", "OpVolumeTest");
109
110 // Get pointer to field entities
111 auto field_ents = data.getFieldEntities();
112 if (auto e_ptr = field_ents[0]) {
113
114 MOFEM_LOG("SYNC", Sev::inform)
115 << "Test " << e_ptr->getEntTypeName() << " " << side;
116
117 // Check if data are cached on entity, and if code is correct, data should
118 // accessible.
119 if (auto ptr = e_ptr->getSharedStoragePtr<MyStorage>()) {
120
121 MOFEM_LOG("SYNC", Sev::verbose)
122 << data.getIndices() << " : " << ptr->globalIndices;
123
124 // Check constancy of data. Stored data are indices, and expected stored
125 // that should be indices, thus difference between should be zero.
126 auto diff = data.getIndices() - ptr->globalIndices;
127 auto dot = inner_prod(diff, diff);
128 if (dot > 0)
129 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "Test failed");
130
131 } else {
132 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
133 "Pointer to stored data is expected to be set");
134 }
135
136 } else {
137 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
138 "Pointer to entity should be set by finite element");
139 }
140
142 }
143};
144
145struct OpVolumeAssemble : public VolOp {
146 OpVolumeAssemble(const std::string &field_name, Vec v)
147 : VolOp(field_name, OPROW), V(v) {}
151 auto nb_dofs = data.getIndices().size();
152 if (nb_dofs) {
153 std::vector<double> v(nb_dofs, 1);
154 CHKERR VecSetValues<EssentialBcStorage>(V, data, &v[0], ADD_VALUES);
155 }
157 }
158
159private:
161};
162
163struct VolRule {
164 int operator()(int, int, int) const { return 1; }
165};
166
167int main(int argc, char *argv[]) {
168
169 // initialize petsc
170 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
171
172 try {
173
174 // Create MoAB database
175 moab::Core moab_core;
176 moab::Interface &moab = moab_core;
177
178 // Create MoFEM database and link it to MoAB
179 MoFEM::Core mofem_core(moab);
180 MoFEM::Interface &m_field = mofem_core;
181
182 // Register DM Manager
183 DMType dm_name = "DMMOFEM";
184 CHKERR DMRegister_MoFEM(dm_name);
185
186 // Simple interface
187 Simple *simple_interface;
188 CHKERR m_field.getInterface(simple_interface);
189 {
190 // get options from command line
191 CHKERR simple_interface->getOptions();
192 // load mesh file
193 CHKERR simple_interface->loadFile();
194 // add fields
195 CHKERR simple_interface->addDomainField("FIELD", H1,
197 // set fields order
198 CHKERR simple_interface->setFieldOrder("FIELD", 4);
199 // setup problem
200 CHKERR simple_interface->setUp();
201
202 // get dm
203 auto dm = simple_interface->getDM();
204 // create elements
205 auto domain_fe = boost::make_shared<VolEle>(m_field);
206 // set integration rule
207 domain_fe->getRuleHook = VolRule();
208
209 auto get_skin = [&]() {
210 Range ents;
211 CHKERR m_field.get_moab().get_entities_by_dimension(0, 3, ents, true);
212 Skinner skin(&m_field.get_moab());
213 Range skin_ents;
214 CHKERR skin.find_skin(0, ents, false, skin_ents);
215 // filter not owned entities, those are not on boundary
216 ParallelComm *pcomm =
217 ParallelComm::get_pcomm(&m_field.get_moab(), MYPCOMM_INDEX);
218 Range proc_skin;
219 CHKERR pcomm->filter_pstatus(skin_ents,
220 PSTATUS_SHARED | PSTATUS_MULTISHARED,
221 PSTATUS_NOT, -1, &proc_skin);
222 Range adj;
223 for (auto d : {0, 1, 2})
224 CHKERR m_field.get_moab().get_adjacencies(proc_skin, d, false, adj,
225 moab::Interface::UNION);
226 // proc_skin.merge(adj);
227 return proc_skin;
228 };
229
230 auto get_mark_skin_dofs = [&](Range &&skin) {
231 auto problem_manager = m_field.getInterface<ProblemsManager>();
232 auto marker_ptr = boost::make_shared<std::vector<unsigned char>>();
233 problem_manager->markDofs(simple_interface->getProblemName(), ROW, skin,
234 *marker_ptr);
235 return marker_ptr;
236 };
237
240 auto mark_dofs = get_mark_skin_dofs(get_skin());
241
242 // set operator to the volume elements
243 domain_fe->getOpPtrVector().push_back(new OpVolumeSet("FIELD"));
244 domain_fe->getOpPtrVector().push_back(new OpVolumeTest("FIELD"));
245
246 domain_fe->getOpPtrVector().push_back(
247 new OpSetBc("FIELD", true, mark_dofs));
248 domain_fe->getOpPtrVector().push_back(new OpVolumeAssemble("FIELD", v));
249 domain_fe->getOpPtrVector().push_back(new OpUnSetBc("FIELD"));
250
251 // make integration in volume (here real calculations starts)
252 CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
253 domain_fe);
255
257 CHKERR VecAssemblyBegin(v);
258 CHKERR VecAssemblyEnd(v);
259 CHKERR VecGhostUpdateBegin(v, ADD_VALUES, SCATTER_REVERSE);
260 CHKERR VecGhostUpdateEnd(v, ADD_VALUES, SCATTER_REVERSE);
261 CHKERR VecGhostUpdateBegin(v, INSERT_VALUES, SCATTER_FORWARD);
262 CHKERR VecGhostUpdateEnd(v, INSERT_VALUES, SCATTER_FORWARD);
263
264 CHKERR DMoFEMMeshToLocalVector(simple_interface->getDM(), v,
265 INSERT_VALUES, SCATTER_REVERSE);
266
267 // CHKERR VecView(v,PETSC_VIEWER_STDOUT_WORLD);
268
269 const MoFEM::Problem *problem_ptr;
270 CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
271
272 double *array;
273 CHKERR VecGetArray(v, &array);
274 for (auto dof : *(problem_ptr->getNumeredRowDofsPtr())) {
275 auto loc_idx = dof->getPetscLocalDofIdx();
276
277 if (dof->getFieldData() != array[loc_idx])
278 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
279 "Data inconsistency");
280
281 if ((*mark_dofs)[loc_idx]) {
282 if (array[loc_idx] != 0)
283 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
284 "Dof assembled, but it should be not");
285 } else {
286 if (array[loc_idx] == 0)
287 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
288 "Dof not assembled, but it should be");
289 }
290 }
291
292 CHKERR VecRestoreArray(v, &array);
293 }
294 }
296
297 // finish work cleaning memory, getting statistics, etc.
299
300 return 0;
301}
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:338
@ 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
@ 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
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:470
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMMoFEM.cpp:373
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:533
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMMoFEM.cpp:1114
MoFEMErrorCode VecSetValues< EssentialBcStorage >(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Set values to vector in operator.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:332
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:277
const double v
phase velocity of light in medium (cm/ns)
const FTensor::Tensor2< T, Dim, Dim > Vec
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
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
UBlasVector< int > VectorInt
Definition: Types.hpp:67
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
CoreTmp< 0 > Core
Definition: Core.hpp:1086
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
constexpr auto field_name
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() 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)
const VectorFieldEntities & getFieldEntities() const
get field entities
const VectorInt & getIndices() const
Get global indices of dofs on entity.
@ OPROW
operator doWork function is executed on FE rows
Set indices on entities on finite element.
keeps basic data about problem
auto & getNumeredRowDofsPtr() const
get access to numeredRowDofsPtr storing DOFs on rows
Problem manager is used to build and partition problems.
MoFEMErrorCode markDofs(const std::string problem_name, RowColData rc, const enum MarkOP op, const Range ents, std::vector< unsigned char > &marker) const
Create vector with marked indices.
Simple interface for fast problem set-up.
Definition: Simple.hpp:26
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:374
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:290
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:806
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:304
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:589
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:749
const std::string getProblemName() const
Get the Problem Name.
Definition: Simple.hpp:334
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:313
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
MyStorage(VectorInt &data)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpVolumeAssemble(const std::string &field_name, Vec v)
Operator set cache stored data, in this example, global indices, but it can be any structure.
boost::shared_ptr< MyStorage > SharedVecInt
OpVolumeSet(const std::string &field_name)
static std::vector< SharedVecInt > entsIndices
This is global static storage.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Test if cached data can be accessed, and check consistency of data.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpVolumeTest(const std::string &field_name)
Set integration rule.
int operator()(int, int, int) const
int main(int argc, char *argv[])
static char help[]