|
From: Laryssa A. <la....@gm...> - 2023-01-06 22:18:14
|
Hello,
I've been having trouble evaluating a cell-centered field (stored in a
ExplicitSystem) at a point using *Number libMesh::System::point_value (
unsigned int var, const Point & p, const Elem * e). *More specifically, I
want to evaluate the field at the centroid of the elements. For most
points/elements, *libMesh::System::point_value* returns zero, but for some
elements it returns the expected non-zero value. Please see the toy example
below.
Does anyone have suggestions on how to deal with this?
Libmesh version: 1.5.1
Thanks,
Laryssa
This is a toy example:
#include "libmesh/exodusII_io.h"
#include "libmesh/mesh.h"
#include "libmesh/elem.h"
#include "libmesh/mesh_base.h"
#include "libmesh/explicit_system.h"
#include "libmesh/equation_systems.h"
#include "libmesh/dof_map.h"
int main(int argc, char ** argv)
{
using namespace libMesh;
LibMeshInit init(argc, argv, MPI_COMM_WORLD);
Mesh mesh_fiber(init.comm());
ExodusII_IO importer(mesh_fiber);
importer.read("file.e");
mesh_fiber.prepare_for_use();
libMesh::EquationSystems es_fiber(mesh_fiber);
typedef libMesh::ExplicitSystem FiberSystem;
FiberSystem& fiber_system = es_fiber.add_system<FiberSystem>("fibers");
fiber_system.add_variable("fibersX", libMesh::CONSTANT, libMesh
::MONOMIAL);
fiber_system.add_variable("fibersY", libMesh::CONSTANT, libMesh
::MONOMIAL);
fiber_system.add_variable("fibersZ", libMesh::CONSTANT, libMesh
::MONOMIAL);
fiber_system.init();
unsigned int step=1;
//Read fiber data into mesh_fiber
importer.copy_elemental_solution(fiber_system, "fibersX", "fibersX",
step);
importer.copy_elemental_solution(fiber_system, "fibersY", "fibersY",
step);
importer.copy_elemental_solution(fiber_system, "fibersZ", "fibersZ",
step);
fiber_system.print_info();
libMesh::MeshBase::const_element_iterator el = mesh_fiber.
active_local_elements_begin();
const libMesh::MeshBase::const_element_iterator end_el = mesh_fiber.
active_local_elements_end();
std::vector<libMesh::Number> ff(3);
for (; el != end_el; ++el)
{
Elem * elem= *el;
libMesh::Point centroid = elem->centroid();
ff[0] = fiber_system.point_value(fiber_system.variable_number(
"fibersX"), centroid, *elem);
ff[1] = fiber_system.point_value(fiber_system.variable_number(
"fibersY"), centroid, *elem);
ff[2] = fiber_system.point_value(fiber_system.variable_number(
"fibersZ"), centroid, *elem);
std::cout << "ff = (" << ff[0] << ", " << ff[1] << ", " << ff[2] <<")
, norm=" << std::sqrt(ff[0]*ff[0]+ff[1]*ff[1]+ff[2]*ff[2]) << std::endl;
}
return 0;
}
|