GridFunction, ProjectCoefficient and GetValue about pymfem HOT 5 OPEN

Without reading ProjectCoefficient, I would guess as follows.
This is probably because H1 element has DoF on mesh node and the projection is done element-wise.
Suppose you have two elements like this
X----0----X----1----X
The left element has attribute=0 and the right one has attribute=1. H1 element DoF is located on node, indicated
as X. The DoFs for the middle X will be determined which element is processed last. If ProjectCoefficient processes
the left one first, and then, the right one, this middle DoF will be 1. If it processes the right one first, then, it will be 0.

If you use L2 element with order=0, you can avoid this issue.

from pymfem.

Hi !

Edit: H1_FECollection has default Gauss-Lobatto rules... which actually makes sense... But when I construct my GridFunctions, I don't seem to be able to evaluate at different quadrature points...

fec = mfem.H1_FECollection(self._order, dim)  # this has Gauss-Lobatto points
scalar_fes = mfem.FiniteElementSpace(mesh, fec)  # This has default values. It should have Lobatto
scalar_fes = mfem.FiniteElementSpace(mesh, fec, mfem.BasisType.GaussLegendre)  # This seems to be ignored.  I get the same thing

thanks for the reply.
Yes, your sketch is what I had in mind. That's why I got confused when imposing Gauss-Legendre rules (although I think they are the default). If the ProjectCoefficient is the operation

$$ x_{node} = \sum_j u_j \cdot \phi_j(node) $$

then the X (interface) should never be evaluated (GL nodes) and x_node should have a unique value (1 or 2). I think I am misunderstanding the projection...

If you use L2 element with order=0, you can avoid this issue.

Actually the binary values (1, 2) were for illustration. I have a density function evaluated at the integration points so constant elements won't work. I also want H1 to avoid jumps at the interface.

All good. If what is happening is that the interfaces have values 1 or 2 depending on which element we are pointing from, I am OK with it. I was just worried that I couldn't reproduce the results because I was doing the wrong thing.

Thanks for the ideas!

from pymfem.

Last question on this: I have tried without success to run the displacement problem with a H1_collection space but then when outputting the solution, evaluate at internal nodes so when I visualize the solution with paraview, each material has the appropriate density value. My problem is that when I do this:

# gf is my gridfunction

        fes = gf.FESpace()
        mesh = fes.GetMesh()
        # nodes = mesh.GetNodes()
        # nodes += gf

        for i in range(fes.GetNE()):
            fe = fes.GetFE(i)
            fdof = fe.GetDof()
            logging.debug(f' TOTAL DOFS {fdof}')
            T = fes.GetElementTransformation(i)
            ir = fe.GetNodes()
            irG = mfem.IntegrationRule(fdof)
            logging.debug(f' INTEGRATION RULE {irG.GetNPoints()}')
            for j in range(fdof):
                ip = ir.IntPoint(j)
                ipG = irG.IntPoint(j)
                logging.debug(f" ORIGINAL {ip.x, ip.y, ip.z, ip.weight}  NOW {ipG.x, ipG.y, ipG.z, ipG.weight}")
                T.SetIntPoint(ip)

the points in the integration rule irG are all set to zero

DEBUG:root: INTEGRATION RULE 4
DEBUG:root: ORIGINAL (0.0, 0.0, 0.0, 0.0)  NOW (0.0, 0.0, 0.0, 0.0)
DEBUG:root: ORIGINAL (1.0, 0.0, 0.0, 0.0)  NOW (0.0, 0.0, 0.0, 0.0)
DEBUG:root: ORIGINAL (0.0, 1.0, 0.0, 0.0)  NOW (0.0, 0.0, 0.0, 0.0)
DEBUG:root: ORIGINAL (0.0, 0.0, 1.0, 0.0)  NOW (0.0, 0.0, 0.0, 0.0)

I know I can set the point ipG = Set3(args) but I am sure there is a direct way to assign it to be Legendre nodes (or whatever). Please, can you point me to the right answer?

Many thanks !

from pymfem.

mfem.IntegrationRule is an array of integration point, and thus, mfem.IntegrationRule(4) will return 4 integration points as array. Perhaps, what you want to do is to use mfem.IntegrationRules to construct a set of rules. Then, ask the rules for an integration rule for specific geometry and order?

>>> rules =  mfem.IntegrationRules(0, mfem.Quadrature1D.GaussLegendre)
>>> rule = rules.Get(1, 2)  ### GeometryType=1, Integration Order=2

Hope this helps.

from pymfem.

Thank you !
Yes, this helps and I can now evaluate at the new points.

only that if I want to write to vtk, I need to update the mesh nodes which are still associated with the GaussLobatto rules...

I need to rethink the approach. I am a bit confused on how to get both things: on one hand, I want the displacement to be C1 (like ex2) even at the interfaces where there is a change of materials, so we use a FE_H1 collection. On the other hand, the young modulus associated to the grid function should not be C1 since it is by definition, a different material. A jump at the change of material interface is expected. But the displacement solution uses the underlying Lame coefficients which are associated to the Young modulus so they have to be in the same space... I am not sure how bonded contacts work anymore.

Thank you for your help.

from pymfem.