147 {
149
150#ifndef NDEBUG
153 << "Refinement for hexes is not implemented";
154#endif
155
156 moab::Core core_ref;
157 moab::Interface &moab_ref = core_ref;
158
159 auto create_reference_element = [&moab_ref]() {
161 constexpr double base_coords[] = {0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0,
162 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1};
164 for (int nn = 0; nn < 8; nn++)
165 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
167 CHKERR moab_ref.create_element(MBHEX, nodes, 8, hex);
169 };
170
171 auto add_ho_nodes = [&]() {
174 CHKERR moab_ref.get_entities_by_type(0, MBHEX, hexes,
true);
176 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
177 CHKERR moab_ref.add_entities(meshset, hexes);
178 CHKERR moab_ref.convert_entities(meshset,
true,
true,
true);
179 CHKERR moab_ref.delete_entities(&meshset, 1);
181 };
182
183 auto set_gauss_pts = [&](std::map<EntityHandle, int> &little_map) {
186 CHKERR moab_ref.get_entities_by_type(0, MBHEX, hexes,
true);
188 CHKERR moab_ref.get_connectivity(hexes, hexes_nodes,
false);
190 gauss_pts.resize(hexes_nodes.size(), 4, false);
191 size_t gg = 0;
192 for (auto node : hexes_nodes) {
193 CHKERR moab_ref.get_coords(&node, 1, &gauss_pts(gg, 0));
194 little_map[node] = gg;
195 ++gg;
196 }
197 gauss_pts = trans(gauss_pts);
199 };
200
201 auto set_ref_hexes = [&](std::map<EntityHandle, int> &little_map) {
204 CHKERR moab_ref.get_entities_by_type(0, MBHEX, hexes,
true);
205 size_t hh = 0;
207 for (auto hex : hexes) {
209 int num_nodes;
210 CHKERR moab_ref.get_connectivity(hex, conn, num_nodes,
false);
211 if (ref_hexes.size2() != num_nodes) {
212 ref_hexes.resize(hexes.size(), num_nodes);
213 }
214 for (int nn = 0; nn != num_nodes; ++nn) {
215 ref_hexes(hh, nn) = little_map[conn[nn]];
216 }
217 ++hh;
218 }
220 };
221
222 auto set_shape_functions = [&]() {
226 const auto nb_gauss_pts = gauss_pts.size2();
227 shape_functions.resize(nb_gauss_pts, 8);
228 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
229 const double ksi = gauss_pts(0, gg);
230 const double zeta = gauss_pts(1, gg);
231 const double eta = gauss_pts(2, gg);
240 }
242 };
243
247
248 CHKERR create_reference_element();
251 std::map<EntityHandle, int> little_map;
252 CHKERR set_gauss_pts(little_map);
253 CHKERR set_ref_hexes(little_map);
254 CHKERR set_shape_functions();
255
257}
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MOFEM_LOG(channel, severity)
Log.
double zeta
Viscous hardening.
std::vector< ublas::matrix< int > > levelRef
std::vector< MatrixDouble > levelShapeFunctions
std::vector< MatrixDouble > levelGaussPtsOnRefMesh