Evaluate field at artbitray position.
std::array<double, 3> point = {0, 0, 0};
std::array<double, 6> eval_points = {-1., -1., -1., 1., 1., 1. };
using VolEle = VolumeElementForcesAndSourcesCore;
auto data_ptr =
getData<VolEle>(point.data(), point.size/3);
if(auto vol_ele = data_ptr->feMethod.lock()) {
vol_ele->getOpPtrVector().push_back(
new MyOp());
auto cache_ptr = boost::make_shared<CacheTuple>();
CHKERR m_field.cache_problem_entities(prb_ptr->getName(), cache_ptr);
"FINITE_ELEMENT_NAME",
data_ptr, m_field.get_comm_rank(),
m_field.get_comm_rank(), cache_ptr);
}
166 CoreInterface &m_field =
cOre;
169 std::vector<EntityHandle> leafs_out;
170 CHKERR data_ptr->treePtr->distance_search(point, distance, leafs_out);
172 for (
auto lit : leafs_out)
173 CHKERR m_field.get_moab().get_entities_by_dimension(lit,
D, tree_ents,
177 MOFEM_LOG(
"SELF", Sev::noisy) <<
"tree entities: " << tree_ents;
179 data_ptr->evalPointEntityHandle.resize(data_ptr->nbEvalPoints);
180 std::fill(data_ptr->evalPointEntityHandle.begin(),
181 data_ptr->evalPointEntityHandle.end(), 0);
183 std::vector<double> local_coords;
184 std::vector<double> shape;
186 auto nb_eval_points = data_ptr->nbEvalPoints;
188 for (
auto tet : tree_ents) {
191 CHKERR m_field.get_moab().get_connectivity(tet, conn, num_nodes,
true);
193 if constexpr (
D == 3) {
195 local_coords.resize(3 * nb_eval_points);
196 shape.resize(4 * nb_eval_points);
198 std::array<double, 12> coords;
199 CHKERR m_field.get_moab().get_coords(conn, num_nodes, coords.data());
202 coords.data(), &data_ptr->evalPoints[0], nb_eval_points,
204 CHKERR Tools::shapeFunMBTET<3>(&shape[0], &local_coords[0],
205 &local_coords[1], &local_coords[2],
210 &shape[0], &shape[1], &shape[2], &shape[3]};
212 &data_ptr->shapeFunctions(0, 0), &data_ptr->shapeFunctions(0, 1),
213 &data_ptr->shapeFunctions(0, 2), &data_ptr->shapeFunctions(0, 3)};
217 &local_coords[0], &local_coords[1], &local_coords[2]};
219 &(data_ptr->localCoords(0, 0)), &(data_ptr->localCoords(0, 1)),
220 &(data_ptr->localCoords(0, 2))};
222 for (
int n = 0;
n != nb_eval_points; ++
n) {
224 const double eps = data_ptr->eps;
225 if (t_shape(0) >= 0 -
eps && t_shape(0) <= 1 +
eps &&
227 t_shape(1) >= 0 -
eps && t_shape(1) <= 1 +
eps &&
229 t_shape(2) >= 0 -
eps && t_shape(2) <= 1 +
eps &&
231 t_shape(3) >= 0 -
eps && t_shape(3) <= 1 +
eps) {
233 data_ptr->evalPointEntityHandle[
n] = tet;
234 t_shape_data(i4) = t_shape(i4);
235 t_local_data(j3) = t_local(j3);
245 if constexpr (
D == 2) {
247 local_coords.resize(2 * nb_eval_points);
248 shape.resize(3 * nb_eval_points);
250 std::array<double, 9> coords;
251 CHKERR m_field.get_moab().get_coords(conn, num_nodes, coords.data());
254 coords.data(), &data_ptr->evalPoints[0], nb_eval_points,
256 CHKERR Tools::shapeFunMBTRI<2>(&shape[0], &local_coords[0],
257 &local_coords[1], nb_eval_points);
261 &shape[0], &shape[1], &shape[2]};
263 &data_ptr->shapeFunctions(0, 0), &data_ptr->shapeFunctions(0, 1),
264 &data_ptr->shapeFunctions(0, 2)};
268 &local_coords[0], &local_coords[1]};
269 data_ptr->localCoords.resize(nb_eval_points, 2);
271 &(data_ptr->localCoords(0, 0)), &(data_ptr->localCoords(0, 1))};
273 for (
int n = 0;
n != nb_eval_points; ++
n) {
275 const double eps = data_ptr->eps;
276 if (t_shape(0) >= 0 -
eps && t_shape(0) <= 1 +
eps &&
278 t_shape(1) >= 0 -
eps && t_shape(1) <= 1 +
eps &&
280 t_shape(2) >= 0 -
eps && t_shape(2) <= 1 +
eps) {
282 data_ptr->evalPointEntityHandle[
n] = tet;
283 t_shape_data(i3) = t_shape(i3);
284 t_local_data(j2) = t_local(j2);
295 const Problem *prb_ptr;
296 CHKERR m_field.get_problem(problem, &prb_ptr);
297 boost::shared_ptr<NumeredEntFiniteElement_multiIndex> numered_fes(
301 in_tets.insert_list(data_ptr->evalPointEntityHandle.begin(),
302 data_ptr->evalPointEntityHandle.end());
303 in_tets = in_tets.subset_by_dimension(
D);
306 MOFEM_LOG(
"SELF", Sev::noisy) <<
"in tets: " << in_tets << endl;
309 m_field.get_finite_elements()->get<FiniteElement_name_mi_tag>().find(
312 m_field.get_finite_elements()->get<FiniteElement_name_mi_tag>().end()) {
314 for (
auto peit = in_tets.pair_begin(); peit != in_tets.pair_end(); ++peit) {
317 prb_ptr->numeredFiniteElementsPtr->get<Unique_mi_tag>().lower_bound(
319 peit->first, (*fe_miit)->getFEUId()));
321 prb_ptr->numeredFiniteElementsPtr->get<Unique_mi_tag>().upper_bound(
323 peit->second, (*fe_miit)->getFEUId()));
324 numered_fes->insert(lo, hi);
327 std::cerr <<
"numered elements:" << std::endl;
328 for (; lo != hi; ++lo)
330 std::cerr << **lo << endl;
333 std::cerr << std::endl;
335 if (
auto fe_ptr = data_ptr->feMethodPtr.lock()) {
338 "Number elements %d to evaluate at proc %d",
339 numered_fes->size(), m_field.get_comm_rank());
343 cache_ptr = boost::make_shared<CacheTuple>();
344 CHKERR m_field.cache_problem_entities(prb_ptr->getName(), cache_ptr);
346 MOFEM_LOG(
"FieldEvaluatorSync", Sev::noisy)
347 <<
"If you call that function many times in the loop consider to "
349 "cache_ptr outside of the loop. Otherwise code can be slow.";
352 CHKERR m_field.loop_finite_elements(prb_ptr, finite_element, *fe_ptr,
353 lower_rank, upper_rank, numered_fes,
354 bh, cache_ptr, verb);
358 "Pointer to element does not exists");