Example of calculating/plot stress or strain on ex2.py/ex17.py?

ddkn commented 1 year ago

Hello,

TL;DR Questions are at the bottom.

I was wondering if anyone had an example or direction for calculating the stress/strain coefficients for ex2.py. I was unsure how to go through generating a new solution (stress or strain) to plot the in GLVis.

I was initially confused and was perusing the C++ mfem issues and came across this, where they define a mfem.CoefficientVector class which has some overloading of methods to calculate the coefficients. I also in turn saw ex17.py for PyMFEM. Just for context here is the GitHub mfem issue C++ a little cleaned up.

/*
 * NOTE according to the issue it should be
 * class StressCoefficient : public CoefficientVector
 * /
class StressCoefficient : public Coefficient
{
private:
   GridFunction &u;
   Vector &s;
   Coefficient &lambda, &mu;
   DenseMatrix eps, sigma;
   int size;
   int i;

public:
   StressCoefficient (GridFunction &_u, Vector &_s, Coefficient &_lambda, Coefficient &_mu, int _size)
      : u(_u), s(_s), lambda(_lambda), mu(_mu), size(_size)
      {
         i=0;
      }
   // void SetComponent(int i, int j) { si = i; sj = j; }
   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
   virtual ~StressCoefficient() { }
};

double StressCoefficient::Eval(ElementTransformation &T,
                               const IntegrationPoint &ip)
{
   u.GetVectorGradient(T, eps);  // eps = grad(u)
   eps.Symmetrize();             // eps = (1/2)*(grad(u) + grad(u)^t)
   double l = lambda.Eval(T, ip);
   double m = mu.Eval(T, ip);
   sigma.Diag(1*eps.Trace(), eps.Size()); // sigma = lambda*trace(eps)*I
   sigma.Add(2*m, eps);          // sigma += 2*mu*eps

   if (eps.Size()== 2) // for 2d
   {
      s[size*0+i]=sigma(0,0); s[size*1+i]=sigma(0,1);
      s[size*2+i]=sigma(1,0); s[size*3+i]=sigma(1,1);
   }
   else if(eps.Size() == 3) // for 3d
   {
      s[size*0+i]=sigma(0,0); s[size*1+i]=sigma(0,1); s[size*2+i]=sigma(0,2);
      s[size*3+i]=sigma(1,0); s[size*4+i]=sigma(1,1); s[size*5+i]=sigma(1,2);
      s[size*6+i]=sigma(2,0); s[size*7+i]=sigma(2,1); s[size*8+i]=sigma(2,2);
   }
   i=i+1;
   return sigma(0,0); // return sig xx to first component of grid function 
}

When I stumbled across ex17.py and found this,

class StressCoefficient(mfem.PyCoefficientBase):
    def __init__(self, lambda_, mu_, si=0, sj=0):
        super(StressCoefficient, self).__init__(0)
        self.lam = lambda_   # coefficient
        self.mu = mu_       # coefficient
        self.si = si
        self.sj = sj     # component
        self.u = None   # displacement GridFunction
        self.grad = mfem.DenseMatrix()

    def SetComponent(self, i, j):
        self.si = i
        self.sj = j

    def SetDisplacement(self, u):
        self.u = u

    def Eval(self, T, ip):
        si, sj = self.si, self.sj
        L = self.lam.Eval(T, ip)
        M = self.mu.Eval(T, ip)
        self.u.GetVectorGradient(T, self.grad)
        if (self.si == self.sj):
            div_u = self.grad.Trace()
            return L * div_u + 2 * M * self.grad[si, si]
        else:
            return M * (self.grad[si, sj] + self.grad[sj, si])

ex17 does have some differences from ex2,

Discontinuous Galerkin (DG) form du/dt + div(sigma(u)) = 0
- where du/dt = 0 which allows a valid solution to use for -div(sigma(u)) = 0
Uses the Gauss-Lobatto basis type
Addition of the DGElasticityIntegrator for InteriorFace and BoundaryFace
A is defined as a SparseMatrix instead of OperatorPtr

Questions

My questions are as follows,

Can you use the definition of $\nabla\cdot u = Trace(\nabla u)$ in example 2 such that

self.u.GetVectorGradient(T, self.grad)
div_u = self.grad.Trace()

Or is this only allow because ex17 is using DG, Gauss-Lobatto basis and a SparseMatrix?

The coefficents class operates on elements (where I see $s_i, s_j$). Is it not possible to use numpy / scipy to do operations? For example, can you do,
```
div_u = numpy.trace(numpy.gradient(u))
```
or equivalent, without using the coefficients class, or is it more work in the end?

I assume it breaks the compatibility with exporting for glvis etc... Or can you do numpy/scipy math and convert it back into equivalent components?

Can you export the projections and plot in matplotlib? or something that has more control for figure labeling and generation?
If you want to isolate a point or area in the mesh to extract values of, how would one go about doing that. Does a boundary need to exist to interact with it? Lets say in beam-tet.mesh you want to get the middle of material 2 on the top part, some area or set of nodes

Figure 1: Boundary example (ex2, ex17) from [https://mfem.org/examples/]

Aside

I may be a bit rusty on my math (a few years has gone by, haha) but here

else:
            return M * (self.grad[si, sj] + self.grad[sj, si])

The $\lambda \nabla\cdot u$ goes to zero when $s_i \neq s_j$. This is because the since the $Trace(\nabla u)$ only has non-zero elements along the trace, we can safely ignore it here?

Sorry for the long post, but I appreciate any guidance!

Cheers,

ddkn commented 1 year ago

In the case for ex17, how would you generate solutions for the stress (*.gf) files given you have the coefficients calculated?

Would it be like this?

stress_c.SetComponent(si, sj)
stress.ProjectCoefficient(stress_c)
stress.Save("filename.gf", 8)

edit: This works! And I wasn't aware you can change the coefficients and change the save said solutions... neat!

ddkn commented 1 year ago

Hm, seems as though I answered for myself (1) and the Aside. This makes sense since in this case u is a unitary symmetric tensor.

For (1) yes you can just use the same solution for plotting the stress on the default basis with the regular ElasticityIntegrator. Is there a problem doing it this way? or issues where I should focus on the Gauss-Lobatto + DGElasticityIntegrator?

~~A new question~~

~~5. With ex17.py, can you visualize the stress on the deformed mesh in GLVis? Or is this bad form? As of now it shows the the stress on the undeformed mesh.~~

Deformed Config

Stress(z, z)

_edit: Nevermind with (5) I realize I can send_solution(mesh, x) then send_solution(mesh, stress)_

ddkn commented 1 year ago

If you want to isolate a point or area in the mesh to extract values of, how would one go about doing that. Does a boundary need to exist to interact with it? Lets say in beam-tet.mesh you want to get the middle of material 2 on the top part, some area or set of nodes

I realize this is pretty easy now, you can use mfem.Mesh.FindPoints in combination with GridFunction.GetValue, for example I have a cylinder rod, where the geometry is the radius $r$ and height $z$ is defined as $r=0.005$, $z=0.5$, you can effectively do this once you have calculated the GridFunction x (or displacement vector) with a known mesh mesh.

xyz = [0.005, 0.0, 0.025]
# Needs a list of coordinates
counts, elem_ids, ips = mesh.FindPoints([xyz])
# Dimensions offset starts at 1
dim = {'x': 1, 'y': 2, 'z': 3}

# Index of elem_ids and ips is for how many "points" you have, in this case 1 so index 0
components = []
components.append(x.GetValue(elem_ids[0], ips[0], dim['x']))
components.append(x.GetValue(elem_ids[0], ips[0], dim['y']))
components.append(x.GetValue(elem_ids[0], ips[0], dim['z']))

magnitude = np.sqrt(sum([i ** 2 for i in components]))

print(components, magnitude)
([-3.1248221669547106e-07, -1.2699223699573933e-09, 4.6536795633482945e-06],
 4.664159112506703e-06)

Which I was able to verify using ParaView

Last question

The only question I have is by using the mfem.PyCoefficientBase you can get the projection of the stress as in ex17, but how does one plot the magnitude such that you can cycle through the projections in GLVis?

I think with that last question I can close the topic, sorry for the noise.

sshiraiwa commented 1 year ago

@ddkn Hi. looks like you are figuring out how to do all by yourself. Great ! Regarding the last point, I am wondering if it is okay to make a new coefficient which returns an absolute value?

ddkn commented 1 year ago

Hi @sshiraiwa, thanks for reaching out!

Hm, I suppose that could work. However, it would be nice to have the gridfunction like the displacement, where in GLVis you can hit F and swap through v_x, v_y, v_z.

I was reading over in mfem/mfem on this issue where they use the VectorCoefficient. I am unsure if this is possible to use this with VectorPyCoefficient, but here is what I have so far?

class StrainVectorCoefficient(mfem.VectorPyCoefficient):
    def __init__(
            self, 
            lambda_coef, 
            mu_coef,):
        # XXX: __init__ cannot be zero, but any number above works
        super(StrainVectorCoefficient, self).__init__(1)
        self.lam = lambda_coef   # coefficient
        self.mu = mu_coef       # coefficient
        self.u = None   # displacement GridFunction
        self.grad = mfem.DenseMatrix()
        self.i: int = 0

    def SetDisplacement(self, u):
        self.u = u

    def Eval(self, v: mfem.Vector, T, ip):
        self.u.GetVectorGradient(T, self.grad)
        self.grad.Symmetrize()

        size = self.grad.Size()
        if size == 2:
            v[0] = self.grad[0, 0]
            v[1] = self.grad[1, 1]
            v[2] = self.grad[0, 1]
        elif size == 3:
            v[0] = self.grad[0, 0]
            v[1] = self.grad[1, 1]
            v[2] = self.grad[2, 2]
            v[3] = self.grad[0, 1]
            v[4] = self.grad[0, 2]
            v[5] = self.grad[1, 2]

scalar_space = mfem.FiniteElementSpace(mesh, fec)
strain = mfem.GridFunction(scalar_space)
strain_coef = StrainVectorCoefficient(lamb_coef, mu_coef)
strain_coef.SetDisplacement(x)

mesh.GetNodes().Assign(reference_nodes)

strain.ProjectCoefficient(strain_coef)
get_xyz(0.005, 0.0, 0.025, mesh, strain)

# ((x, y, z), magnitude)
([-6.132483942100051e-05, -6.132483942100051e-05, -6.132483942100051e-05],
 0.00010621773764317564)

Now if I compare to changing projection along $(x, x)$, $(y, y)$, $(z, z)$ you can see we get one set of values but not the others. I think I am only getting $x$, not $y$ and $z$.

for i in range(3):
    scalar_space = mfem.FiniteElementSpace(mesh, fec)
    strain = mfem.GridFunction(scalar_space)
    strain_coef = StrainCoefficient(lamb_coef, mu_coef)
    strain_coef.SetDisplacement(x)
    strain_coef.SetComponent(i, i)
    strain.ProjectCoefficient(strain_coef)
    print(i, get_xyz(0.005, 0.0, 0.025, mesh, strain))
# ((x, y, z), magnitude)
0
([-6.132483942100051e-05, -6.132483942100051e-05, -6.132483942100051e-05], 0.00010621773764317564)

1
([-6.130708396019692e-05, -6.130708396019692e-05, -6.130708396019692e-05], 0.00010618698428295205)

2
([0.0001886661919339077, 0.0001886661919339077, 0.0001886661919339077], 0.0003267794301000696)

I would think I would get ([-6.132483942100051e-05, -6.130708396019692e-05, 0.0003267794301000696], 0.0003380888794536733).

ddkn commented 1 year ago

Is there a way to combine the solutions together? To rebuild a gridfunction that has all 3 components?

ddkn commented 1 year ago

Another thing is the VectorCoefficient is projecting on $(x, x)$, where when I plot Strain $(x, x)$, they are the same, if I plot *Strain $(z, z)$ it is different (obviously, haha). But how do I project with the VectorCoefficient?

VectorCoefficient

Coefficient (x, x)

Coefficient (z, z)

sshiraiwa commented 1 year ago

To combine several GridFunctions to one, you can use vdim in GridFunction. From Python, you can do something like this...

  mesh = mfem.mesh.Mesh().MakeCartesian2D(10, 10, 3)
  fec1 = mfem.H1_FECollection(1, 2)   # order 1, dim = 2
  fes1 = mfem.FiniteElementSpace(mesh, fec1)       # 1 component
  fes2 = mfem.FiniteElementSpace(mesh, fec1, 2)   # 2 component
  g1x = mfem.GridFunction(fes1)
  g1y = mfem.GridFunction(fes1)
  g2 = mfem.GridFunction(fes2)
  g1x.Size()    # print 121...
  g3.Size()    # print 242...

  #Let' say you have a data on g1x and g1y
  g1x.GetDataArray()[:] = 1         # set g1x using GetDataArray which is a numpy array pointing the same memory
  g1y.Assign(np.arange(121.))   # you can also do this.

  # then this is to set x and y component to g2
  g2.Assign(np.hstack([g1x.GetDataArray(), g1y.GetDataArray()])

 # to check it....  
  mesh.Save('test.mesh')
  g2.Save('test.gf')

  glvis -m test.mesh -g test.gf -gc 0   ## to plot 0 component

ddkn commented 1 year ago

Thanks @sshiraiwa! That is neat! In here,

#Let' say you have a data on g1x and g1y
g1x.GetDataArray()[:] = 1         # set g1x using GetDataArray which is a numpy array pointing the same memory
g1y.Assign(np.arange(121.))   # you can also do this.

Are you just assigning data to these gridfunctions? Just making sure. In the case I am curious about I would really want to change

fes2 = mfem.FiniteElementSpace(mesh, fec1, 2)   # 2 component

fes3 = mfem.FiniteElementSpace(mesh, fec1, 3)

Where I would just append the z-component to the np.hstack. I could in principle do operations like square, sum and squareroot to get the magnitude on a 4th component gridfunction no?

In a slightly related note, do you know how to change which vector the VectorPyCoefficient class points along as my above issue with it defaulting to $(x, x)$ projection? How would one change to say a $(x, y)$ projection?

sshiraiwa commented 1 year ago

Hi @ddkn, Assign calls copy assignment operator in MFEM C++ (= operator). Thus, it copies the data. In other words, at the step of calling g2.Assign above, (I think) it creates a temporary array. If you want to avoid this copying, you need to get the access to the memory hold by g2 using g2.GetDataArray(). This returns a numpy array pointing to the same memory area. I am not sure about VectorPyCoefficient. Isn't fec H1 element with vdim=1? If so, I suppose that strain in your code is scalar, and it can not hold three values for x, y, and z. Am I missing something here?

ddkn commented 1 year ago

I am not sure about VectorPyCoefficient. Isn't fec H1 element with vdim=1? If so, I suppose that strain in your code is scalar, and it can not hold three values for x, y, and z. Am I missing something here?

Ah I am an idiot! I wasn't paying attention, when I converted from PyCoefficient to VectorPyCofficent. Thanks for catching that! I think I am a bit overwhelmed on the tooling as I am teaching myself FEM.

My code should have changed from,

scalar_space = mfem.FiniteElementSpace(mesh, fec)
strain = mfem.GridFunction(scalar_space)

to

dim = 3 # Which I could pull from when I defined the displacement x
vector_space = mfem.FiniteElementSpace(mesh, fec, dim)
strain = mfem.GridFunction(vector_space)

This change alone fixes all my issues with Last Question. I can now view the $magnitude$, $v_x$, $v_y$, and $v_z$ using F in GLVis!

I guess I didn't realize the mfem.FiniteElementSpace defaulted to 1 dimension (scalar), if this option wasn't passed. Looking at the ouptut from mfem.FiniteElementSpace??, the __init__ comment showed it was overloaded based on parameters. This makes getting stress trivial, now that I understand how this works better. Now on to more fun!

I think this answers all of my questions, and can mark this closed. Thanks again @sshiraiwa, your help was invaluable!

mfem / PyMFEM