-
Notifications
You must be signed in to change notification settings - Fork 62
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Example of calculating/plot stress or strain on ex2.py/ex17.py? #184
Comments
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?
edit: This works! And I wasn't aware you can change the coefficients and change the save said solutions... neat! |
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
edit: Nevermind with (5) I realize I can |
@ddkn Hi. looks like you are figuring out how to do all by yourself. Great ! |
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 I was reading over in mfem/mfem on this issue where they use the 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 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 |
Is there a way to combine the solutions together? To rebuild a gridfunction that has all 3 components? |
To combine several GridFunctions to one, you can use vdim in GridFunction. From Python, you can do
|
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 In a slightly related note, do you know how to change which vector the |
Hi @ddkn, Assign calls copy assignment operator in MFEM C++ (= operator). Thus, it copies the data. |
Ah I am an idiot! I wasn't paying attention, when I converted from 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 I guess I didn't realize the I think this answers all of my questions, and can mark this closed. Thanks again @sshiraiwa, your help was invaluable! |
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.When I stumbled across ex17.py and found this,
ex17 does have some differences from ex2,
SparseMatrix
instead ofOperatorPtr
Questions
My questions are as follows,
Or is this only allow because ex17 is using DG, Gauss-Lobatto basis and a
SparseMatrix
?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
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,
The text was updated successfully, but these errors were encountered: