GridFunction, ProjectCoefficient and GetValue

[XuliaDS]

Hi:
I am trying to use the ElasticIntegrator for a multi-material problem where some of the materials vary the Young modulus as a function of density. So I need to change the way lame constants are computed. In example 2, they are computed using PWConstantCoefficient:

lamb = mfem.Vector(mesh.attributes.Max())
lamb.Assign(1.0)
lamb[0] = lamb[1] * 50
lambda_func = mfem.PWConstCoefficient(lamb)

My approach was to create a LinearCoefficient class which has information about another field (in this case: density). But I have an inconsistency of values when I evaluate the coefficient class and the grid function. I have made a dumb example where this problem arises.

'''
   MFEM example 2

   See c++ version in the MFEM library for more detail
'''
import os
from os.path import expanduser, join, dirname
import numpy as np

import mfem.ser as mfem

class LinearCoefficient(mfem.PyCoefficientBase):
    """
    This is a local class that uses a grid function to get the projection values
    """

    def __init__(self, vals):
        super(LinearCoefficient, self).__init__(0)
        self._v = [v for v in vals]
        self._gf = None

    def SetGF(self, gf):
        self._gf = gf

    def Eval(self, T, ip):
        if self._gf:
            res1 = self._gf.GetValue(T, ip)
            if res1 != self._v[T.Attribute - 1]:
                print(f' ERROR !! we should have the same value at each material but {res1} not { self._v[T.Attribute - 1]}!!! ')
                raise "LEAVE NOW"
        return self._v[T.Attribute - 1]


def run(order=1, meshfile=''):
    '''
    run ex2
    '''

    mesh = mfem.Mesh(meshfile, 1, 1)
    dim = mesh.Dimension()
    fec = mfem.H1_FECollection(order, dim)
    scalar_fes = mfem.FiniteElementSpace(mesh, fec)

    # we have two materials. Assign value = 1, 2
    vals = [1,2]
    myCoeff = LinearCoefficient(vals)
    gf = mfem.GridFunction(scalar_fes)
    gf.ProjectCoefficient(myCoeff)
    myCoeff.SetGF(gf)
    gf.ProjectCoefficient(myCoeff)

if __name__ == "__main__":
    from mfem.common.arg_parser import ArgParser

    parser = ArgParser(description='Ex2 (Linear elasticity)')
    parser.add_argument('-m', '--mesh',
                        default='beam-tri.mesh',
                        action='store', type=str,
                        help='Mesh file to use.')
    parser.add_argument('-o', '--order',
                        action='store', default=1, type=int,
                        help="Finite element order (polynomial degree) or -1 for isoparametric space.")
    args = parser.parse_args()
    parser.print_options(args)

    order = args.order

    meshfile = expanduser(
        join(dirname(__file__), '..', 'data', args.mesh))

    run(order=order, meshfile=meshfile )

The first time we perform the operation gf.ProjectCoefficient(myCoeff) we are assigning a constant value of 1 or 2 according to the element Attribute. Then I store the gridFunction and call it again, this time checking that values coincide.
However, when I run it I get an error:

    raise "LEAVE NOW"
TypeError: exceptions must derive from BaseException
Options used:
   --mesh  beam-tri.mesh
   --order  1
 ERROR !! we should have the same value at each material but 2.0 not 1!!! 

meaning that the grid function at (T, ip) gave different value than the element value. My first thought was that it might be using a quadrature rule including endpoints (-,1 1) but I forced it to GaussLegendre and I get the same error. What am I missing?? Shouldn-t both produce the same result ?

Many thanks

[sshiraiwa]

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.

[XuliaDS]

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!

[XuliaDS]

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 !

[sshiraiwa]

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.

[XuliaDS]

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.