v0.14.0
Loading...
Searching...
No Matches
field_blas.cpp
Go to the documentation of this file.
1/** \file field_blas.cpp
2 * \example field_blas.cpp
3 * \brief test field blas interface
4 *
5 * \ingroup mofem_field_algebra
6 */
7
8
9
10#include <MoFEM.hpp>
11
12using namespace MoFEM;
13
14static char help[] = "...\n\n";
15
16int main(int argc, char *argv[]) {
17
18 MoFEM::Core::Initialize(&argc, &argv, PETSC_NULL, help);
19
20 try {
21
22 moab::Core mb_instance;
23 moab::Interface &moab = mb_instance;
24 int rank;
25 MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
26
27 // Reade parameters from line command
28 PetscBool flg = PETSC_TRUE;
29 char mesh_file_name[255];
30#if PETSC_VERSION_GE(3, 6, 4)
31 CHKERR PetscOptionsGetString(PETSC_NULL, "", "-my_file", mesh_file_name,
32 255, &flg);
33#else
34 CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-my_file",
35 mesh_file_name, 255, &flg);
36#endif
37 if (flg != PETSC_TRUE) {
38 SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
39 }
40 int order = 1;
41#if PETSC_VERSION_GE(3, 6, 4)
42 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-my_order", &order, PETSC_NULL);
43#else
44 CHKERR PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-my_order", &order,
45 PETSC_NULL);
46#endif
47
48 // Read mesh to MOAB
49 const char *option;
50 option = "";
51 CHKERR moab.load_file(mesh_file_name, 0, option);
52
53 // Create MoFEM database
54 // Note: Is MoFEM::CoreTmp<1> for testing purposes only
55 MoFEM::CoreTmp<-1> core(moab);
56 MoFEM::Interface &m_field = core;
57
58 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
59 0, 3, BitRefLevel().set(0));
60
61 // Create fields, add entities to field and set approximation order
62 CHKERR m_field.add_field("FIELD_A", H1, AINSWORTH_LEGENDRE_BASE, 3,
63 MB_TAG_DENSE);
64 CHKERR m_field.add_field("FIELD_B", H1, AINSWORTH_LEGENDRE_BASE, 3,
65 MB_TAG_DENSE);
66 CHKERR m_field.add_field("FIELD_C", H1, AINSWORTH_LEGENDRE_BASE, 1,
67 MB_TAG_DENSE);
68
69 CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "FIELD_A");
70 CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "FIELD_B");
71 CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "FIELD_C");
72
73 CHKERR m_field.set_field_order(0, MBTET, "FIELD_A", order + 1);
74 CHKERR m_field.set_field_order(0, MBTRI, "FIELD_A", order + 1);
75 CHKERR m_field.set_field_order(0, MBEDGE, "FIELD_A", order + 1);
76 CHKERR m_field.set_field_order(0, MBVERTEX, "FIELD_A", 1);
77
78 CHKERR m_field.set_field_order(0, MBTET, "FIELD_B", order);
79 CHKERR m_field.set_field_order(0, MBTRI, "FIELD_B", order);
80 CHKERR m_field.set_field_order(0, MBEDGE, "FIELD_B", order);
81 CHKERR m_field.set_field_order(0, MBVERTEX, "FIELD_B", 1);
82
83 CHKERR m_field.set_field_order(0, MBTET, "FIELD_C", order);
84 CHKERR m_field.set_field_order(0, MBTRI, "FIELD_C", order);
85 CHKERR m_field.set_field_order(0, MBEDGE, "FIELD_C", order);
86 CHKERR m_field.set_field_order(0, MBVERTEX, "FIELD_C", 1);
87
88 // build field
89 CHKERR m_field.build_fields();
90
91 // get access to field agebra interface
92 auto fb = m_field.getInterface<FieldBlas>();
93
94 auto check_axpy = [&]() {
96 // set value to field
97 CHKERR fb->setField(+1, MBVERTEX, "FIELD_A");
98 CHKERR fb->setField(-2, MBVERTEX, "FIELD_B");
99
100 // FIELD_A = FIELD_A + 0.5 * FIELD_B
101 CHKERR fb->fieldAxpy(+0.5, "FIELD_B", "FIELD_A");
102 // FIELD_B *= -0.5
103 CHKERR fb->fieldScale(-0.5, "FIELD_B");
104
105 auto dofs_ptr = m_field.get_dofs();
106 for (auto dit : *dofs_ptr) {
107
108 auto check = [&](const std::string name, const double expected) {
110 if (dit->getName() == name) {
111 cout << name << " " << dit->getFieldData() << " " << expected
112 << endl;
113 if (dit->getFieldData() != expected)
114 SETERRQ2(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
115 "Wrong DOF value 0 != %4.3e for %s", dit->getFieldData(),
116 boost::lexical_cast<std::string>(*dit).c_str());
117 }
118
120 };
121
122 CHKERR check("FIELD_A", 0);
123 CHKERR check("FIELD_B", 1);
124 }
126 };
127
128 auto check_set_vertex = [&]() {
130
131 auto set_distance = [&](VectorAdaptor &&field_data, double *xcoord,
132 double *ycoord, double *zcoord) {
134 FTensor::Index<'i', 3> i;
136 xcoord, ycoord, zcoord);
137 field_data[0] = sqrt(t_coord(i) * t_coord(i));
139 };
140
141 CHKERR fb->setVertexDofs(set_distance, "FIELD_A");
142
143 // Test set veryex
144 struct TestMethod : EntityMethod {
145 TestMethod(moab::Interface &moab) : EntityMethod(), moab(moab) {}
146
147 MoFEMErrorCode preProcess() { return 0; }
150 if (entPtr->getEntType() == MBVERTEX) {
151 EntityHandle v = entPtr->getEnt();
152 VectorDouble3 coords(3);
153 CHKERR moab.get_coords(&v, 1, &coords[0]);
154 if (std::abs(entPtr->getEntFieldData()[0] - norm_2(coords)) >
155 std::numeric_limits<double>::epsilon())
156 SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
157 "Wrong vals %3.4e != %3.4e",
158 entPtr->getEntFieldData()[0], norm_2(coords));
159 }
161 }
162 MoFEMErrorCode postProcess() { return 0; }
163
164 private:
165 moab::Interface &moab;
166 };
167
168 TestMethod method(moab);
169 // checking if all is ok
170 CHKERR m_field.loop_entities("FIELD_A", method);
171
173 };
174
175 auto check_lambda = [&]() {
177
178 Range meshset_ents;
180 m_field, SIDESET | PRESSURESET, it)) {
181 if (it->getMeshsetId() != 1)
182 continue;
183 Range ents, nodes;
184 CHKERR it->getMeshsetIdEntitiesByType(m_field.get_moab(), MBTRI, ents,
185 true);
186 CHKERR m_field.get_moab().get_connectivity(ents, nodes, true);
187 meshset_ents.merge(nodes);
188 }
189 // set value to field
190 CHKERR fb->setField(+1, MBVERTEX, "FIELD_A");
191 CHKERR fb->setField(-1, MBVERTEX, "FIELD_B");
192 CHKERR fb->setField(1.5, MBVERTEX, "FIELD_C");
193
194 auto field_axpy = [&](const double val_y, const double val_x) {
195 return 0.5 * val_x + val_y;
196 };
197 auto field_scale = [&](const double val) { return -0.5 * val; };
198
199 // note that y_ent is first
200 // FIXME: this can be confusing?
201 auto vector_times_scalar_field =
202 [&](boost::shared_ptr<FieldEntity> ent_ptr_y,
203 boost::shared_ptr<FieldEntity> ent_ptr_x) {
205 auto x_data = ent_ptr_x->getEntFieldData();
206 auto y_data = ent_ptr_y->getEntFieldData();
207 // const auto size_x = x_data.size(); // scalar
208 const auto size_y = y_data.size(); // vector
209
210 for (size_t dd = 0; dd != size_y; ++dd)
211 y_data[dd] *= x_data[0];
212
213 // y_data *= x_data[0]; //would work as well
214
216 };
217
218 Range ents; // TODO: create test for subentities
219 // FIELD_A = FIELD_A + 0.5 * FIELD_B
220 CHKERR fb->fieldLambdaOnValues(field_axpy, "FIELD_B", "FIELD_A", true);
221 // CHKERR fb->fieldAxpy(+0.5, "FIELD_B", "FIELD_A");
222
223 // FIELD_B *= -0.5
224 // CHKERR fb->fieldScale(-0.5, "FIELD_B");
225 CHKERR fb->fieldLambdaOnValues(field_scale, "FIELD_B", &meshset_ents);
226
227 // FIELD_B(i) *= FIELD_C
228 CHKERR fb->fieldLambdaOnEntities(vector_times_scalar_field, "FIELD_C",
229 "FIELD_B", true, &meshset_ents);
230
231 auto dofs_ptr = m_field.get_dofs();
232 constexpr double eps = std::numeric_limits<double>::epsilon();
233 for (auto dit : *dofs_ptr) {
234 auto check = [&](const std::string name, const double expected) {
236 if (dit->getName() == name && dit->getEntType() == MBVERTEX) {
237 cout << name << " " << dit->getFieldData() << " " << expected
238 << endl;
239 if (abs(dit->getFieldData() - expected) > eps)
240 SETERRQ3(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
241 "Wrong DOF value %4.3e != %4.3e for %s",
242 dit->getFieldData(), expected,
243 boost::lexical_cast<std::string>(*dit).c_str());
244 }
245
247 };
248
249 CHKERR check("FIELD_A", 0.5);
250
251 if (meshset_ents.find(dit->getEnt()) != meshset_ents.end()) {
252 CHKERR check("FIELD_B", 0.75);
253 CHKERR check("FIELD_C", 1.5);
254 }
255 }
256
258 };
259
260 CHKERR check_axpy();
261 CHKERR check_set_vertex();
262 CHKERR check_lambda();
263 }
265
267 return 0;
268}
int main()
static const double eps
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ H1
continuous field
Definition definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ PRESSURESET
@ SIDESET
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
constexpr int order
static char help[]
virtual const DofEntity_multiIndex * get_dofs() const =0
Get the dofs object.
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_entities(const std::string field_name, EntityMethod &method, Range const *const ents=nullptr, int verb=DEFAULT_VERBOSITY)=0
Loop over field entities.
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
char mesh_file_name[255]
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
VectorBoundedArray< double, 3 > VectorDouble3
Definition Types.hpp:92
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition Types.hpp:115
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
virtual MoFEMErrorCode postProcess()
function is run at the end of loop
virtual MoFEMErrorCode operator()()
function is run for every finite element
virtual MoFEMErrorCode preProcess()
function is run at the beginning of loop
Managing BitRefLevels.
virtual moab::Interface & get_moab()=0
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.
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 structure to exchange data between mofem and User Loop Methods on entities.
EntityMethod()=default
Basic algebra on fields.
Definition FieldBlas.hpp:21
MoFEMErrorCode setField(const double val, const EntityType type, const std::string field_name)
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.