Hi, I would first verify that the System's solution vector has the expected values after calling "copy_elemental_solution". You could do this by just printing out the solution vector if it is a small enough problem. Once you are sure that part is working, I would have to assume that there is some kind of bug with System::point_value() in libmesh-1.5.1. Since that libmesh version is from 2019 ( https://github.com/libMesh/libmesh/wiki/libMesh-release-history), it's unlikely that anyone will be able to spend much time debugging it. If you find that we have the same bug in libmesh master or libmesh-1.7.1, someone might be able to take a closer look.
Also, you may find that you get faster/more varied responses by starting a Discussion on GitHub (https://github.com/libMesh/libmesh/discussions) than by using this old mailing list. -- John On Fri, Jan 6, 2023 at 4:18 PM Laryssa Abdala <la.abd...@gmail.com> wrote: > 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; > } > _______________________________________________ Libmesh-users mailing list Libmesh-users@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/libmesh-users
