Hi Robert:
I’ve been able to get the stiffness matrix under your help, but I still get several puzzles about the linear system components [1] (8):
Is the linear system in the matrix formula: Kx=r ? I exported both K and r to matlab to used the right division to get x: x = r \ K. But this gave me a different result from variables = pb.solve() and x = variables.vec. Analysis showed that pb.solve() gave the appropriate result, what’s the difference between the two approaches?
In my project, the …
[View More]measurements (knowns) are Φ on surface (∂Ω) and the unknown is the source distribution Q in attached ProblemDescription.docx. So I need to separate the vectors Φ and Q (nodal values of Φ(x) and Q(x)) explicitly in the matrix form, which means I need to get the matrix F: (FQ is the discretized form of ∫_Ω▒Q(x)s(x)dx in (2), where Q(x) is the source distribution, s(x) is the test function). How can I get this matrix?
I found that the solver in SfePy was much faster than the right division in matlab, so I’m curious that which function does SfePy adopts to solve the linear system. Besides, my project requires the inversion of a sub-matrix of the stiffness matrix K, does SfePy offer such methods or which function do you recommend to invert such a banded sparse matrix?
[View Less]
Thanks for your patient instructions, and I have solved the former problems under your help. But I’m sorry to have to bother you again.
I could get the stiffness matrix with expected shape now, but the matrix elements are all zeros! I guess it’s most likely to go wrong in the material definition and term selection as well as term implementation steps. I have attached my code snippet together with the mesh used and a description file. Could you please help me check what’s wrong with my code …
[View More]snippet, especially my material definition part, and see if I have chosen the right terms to build equations?
Thanks!
[View Less]
Hi, thank you for the SfePy. I wish to apply SfePy in my research project, but I'm now having a few problems.
I want to get the stiffness matrix for postprocessing by my self. The pb (Problem) object has an attribute namely mtx_a, is it the stiffness matrix I'm looking for? If so, I get some puzzles: in my view, after "pb = Problem('name',equations=eqs)" and "pb.set_bcs", the problem has been already clearly defined. But so far the pb.mtx_a is still a non-type object and I can get the pb.…
[View More]_mtx_a only after "pb.set_solver(nls)", "pb.solve()". I believe the stiffness matrix is determined only by the problem itself and corresponding essential boundary conditions, not including the solvers. So could you please tell me why is that, or where can I find the proper stiffness of the problem?
Looking forward for your answers, and thank you again for this wonderful project!
[View Less]
Got it, thanks a million! But I just got in some new troubles lol. I created a mesh file namely ‘test_mouse.mesh’. I was able to load it, but when I do “domain = FEDomain(‘domain’, mesh)”, a warning occurs: bad element orientation, trying to correct, … corrected. But still I could use this domain to create regions, and resview showed the regions were just as expected. So what does this warning mean, does SfePy requires special vertices indexing? And what happened when “trying to correct“, did …
[View More]it change my previous vertices and/or tetrahedra indexing?
I have attached the mesh file, could please check what is wrong with my mesh? Thanks a lot!
[View Less]
Hi,
thank you for SfePy. I would like to use it for checking an analytic implementation of a linear elastic poro-elasticity problem. However, I am a total beginner if it comes to using FEM methods in general, so this email is a request for help, partially probably also on very basic questions. If there's interest I would be happy to write up my journey as a tutorial for absolute beginners.
The problem that I want to solve is the following:
I want to look at a cuboidal inclusion in a …
[View More]stratified elastic half-space. Therefore, we use the usual momentum balance equation:
0 = grad(sigma) + F
and we use stress-strain relations of the kind
sigma_ij = C_ijkl(e_kl)
We can decompose the strain tensor in two parts: e = e_eff - e_I where e_eff is the effective part and e_I is the part stemming from the inclusion.
the stress tensor due to the inclusion is given by
sigma_ij,I = C_ijkl e_kl,I
We want to consider hydrostatic inclusions, therefore the inclusion stress tensor is diagonal: sigma_ij,I = p * delta_ij.
This strain inclusion produces the same deformation as the external force F = -grad p.
We assume that the inclusion is embedded in a stratified half-space: we have a free surface, and sections parallel to the surface with differing linear elastic isotropic material constants.
Now, with regards to solving this problem with SfePy, I have the following questions:
1. I would like to solve the problem on a hexahedral mesh of rectangular cubes. Creating the vertices, and connections is quite straightforward without the need for a mesh generator. However, it is not obvious to me, how to hand this mesh to SfePy:
I assume that I need to use Mesh.from_data(), furthermore, I assume that my vertices go under coors, and my connectivity matrix goes under conns, and that descs should be set to '3_8', but I am not sure about the other parameters.
2. In connection to that, is there a preferred "Dimensionality” for SfePy? I would usually define my grid in the thousands of m, my elastic constants in Pa (which makes for numbers in the order of 1e9), and which results in deformations in the order of cm. Does this cause numerical issues?
3. I had a look through the predfined equations and linear elasticity examples. It appears, as if I should use the
dw_lin_elastic
term for the force balance, but how do I encode the strain inclusion in the weak formalism?
One option would be to use that my force term is given by F=-grad p. However, if I assume a constant pressure change in my inclusion, this gives a pretty exciting delta function like derivative. Is there
a better way of doing this?
4. In principle, I would like the top boundary of the box as a free surface, and the other ones should essentially "slide": so not being fixed to u=0. Is there a possibility to accomplish this?
5. Which sort of solver would you recommend for this problem?
6. If I would also like to get the stresses and strains, would you recommend to calculate them afterwards from the displacment or is there a better way of getting them?
7. Obviously, the problem is symmetric around the x = 0 and the y = 0 axis. Is it beneficial to exploit this?
Thanks in advance,
Kind regards,
Hanno
P.S. in case you are interested, I am trying to check my impementation of https://doi.org/10.1111/j.1365-246X.2007.03487.x : B. Kvshinov, “Reflectivity method for geomechanical equilibria".
Below is what I have so far (mostly cobbled together from various turorials):
import sys
import numpy as np
import matplotlib.pyplot as plt
from sfepy.mechanics.matcoefs import stiffness_from_youngpoisson
from sfepy.discrete.fem import Mesh, FEDomain
from sfepy.discrete import FieldVariable, Material, Integral, Function, \
Equation, Equations, Problem
from sfepy.discrete.conditions import Conditions, EssentialBC
from sfepy.base.base import IndexedStruct
from sfepy.solvers.ls import ScipyDirect
from sfepy.solvers.nls import Newton
def create_mesh(x0=-4000., y0=-4000., z0=0.,
x1=4000, y1=4000, z1=2000,
dx=100, dy=100, dz=20):
# create the mesh
z_spacing = dz
x_spacing = dx
y_spacing = dy
the_grid = np.mgrid[x0:x1+1:x_spacing,y0:y1+1:y_spacing,z0:z1+1:z_spacing]
xr = the_grid[0].ravel()
yr = the_grid[1].ravel()
zr = the_grid[2].ravel()
vertices = np.vstack((xr, yr, zr)).T
nx, ny, nz = the_grid.shape[1:]
tetrahedra = []
dxs = np.array([ x_spacing, 0., -x_spacing,0., x_spacing,0., -x_spacing])
dys = np.array([0., y_spacing, 0, -y_spacing, 0, y_spacing, 0])
dzs = np.array([0, 0,0,z_spacing,0,0,0])
for i, v in enumerate(vertices):
tet = [i, i+nz*ny, i+nz*ny+nz, i+nz,
i+1, i+nz*ny+1, i+nz*ny+nz+1, i+nz+1]
try:
cube = vertices[tet]
dc = cube[1:] - cube[:-1]
if (np.array_equal(dc[:,0],dxs) & np.array_equal(dc[:,1],dys) &
np.array_equal(dc[:,2],dzs)):
tetrahedra.append(tet)
except IndexError:
continue
tetrahedra = np.asarray(tetrahedra)
return vertices, tetrahedra
descs3d = ['3_8']
x0 = -4000.
x1 = 4000.
dx = 100.
y0 = -4000.
y1 = 4000.
dy = 100.
z0 = 0.
z1 = 2000.
dz = 20.
vertices, tets = create_mesh(x0, y0, z0, x1, y1, z1, dx, dy, dz)
# How do I do this?
# mesh3d = Mesh.from_data('mesh', coors=vertices,
# ngroups=1, conns=tets,
# Define the regions
domain = FEDomain('domain', mesh3d)
# inclusion extent
xr0 = -501.
xr1 = 501.
yr0 = -251.
yr1 = 251.
zr0 = 1000.-1.
zr1 = 1101.
rexp = f'vertices in ((x>{xr0})&(x<{xr1})&(y>{yr0})&(y<{yr1})&(z>{zr0})&(z<{zr1}))'
regions = {'Omega' : 'all',
'free_surface' : ('vertices in (z < 0.01)', 'facet'),
'base_boundary': ('vertices in (z > %f.2)'%(z1 - 0.5), 'facet'),
'Gamma_left': ('vertices in (x < %f.1)'%(x0 + 0.5), 'facet'),
'Gamma_right': ('vertices in (x > %f.1)'%(x1 - 0.5), 'facet'),
'Gamma_front': ('vertices in (y < %f.1)'%(y0+0.5), 'facet'),
'Gamma_back': ('vertices in (y > %f.1)'%(x1-0.5), 'facet'),
'reservoir_layer': (f'vertices in ((z > {zr0}) and (z<{zr1}','cell'),
'reservoir' : (rexp, 'cell'),
'overburden' : ('vertices in (z<{zr0})', 'cell'),
'underburden' : ('vertices in (z>{zr1})', 'cell')}
approx_order = 1
fields = {'displacement' : ('real', 3, 'Omega', approx_order)}
# What do I do about the inclusion here?
materials={'overburden' : ({'D':stiffness_from_youngpoisson(3, 1e9, 0.4)}),
'res_layer' : ({'D':stiffness_from_youngpoisson(3, 2e9, 0.25)}),
'underburden' : ({'D':stiffness_from_youngpoisson(3, 5e9, 0.3)})}
variables = {
'u' : ('unknown field', 'displacement', 0),
'v' : ('test field', 'displacement', 'u')}
# Not sure what to do here.
bcs = {}
integrals = {
'i': 2,
}
##
# Force Balance :
equations = {
'elasticity' : """dw_lin_elastic.i.Omega( solid.D, v, u ) = 0""",}
[View Less]