Example of calculating/plot stress or strain on ex2.py/ex17.py? about pymfem HOT 11 CLOSED

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?

from pymfem.

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!

from pymfem.

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)

from pymfem.

4. 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.

from pymfem.

@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?

from pymfem.

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).

from pymfem.

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

from pymfem.

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)

from pymfem.

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

from pymfem.

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?

from pymfem.

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!

from pymfem.