Re: [Libmesh-users] libMesh::System::point_value returning zero values

[John P. <jwp...@gm...>]

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....@gm...> 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;
> }
>