Hi there,
I followed the interactive example in sfepy tutorial and wrote my custom simulation script for the block domain with body force (t2 term). That works. But for the point loaded block model (t3 term), I can't make it work. The script stopped at "assert_(mat.shape[-1]==virtual.dim)" of terms_point.py with IndexError: tuple index out of range. Searching the web does not find help information for this.
The reason I write the attached custom simulation script is that I failed doing the things below: the "non-custom" script for a point loaded block running with "simple.py" works and gmsh can mesh different geometries. I would like to finish those two things together in python like "os.system('gmsh ........')", and "os.system('simple.py examples/myscript.py") but errors to import config in simple.py got me there since "path" settings for import in simple.py seems a little complicated for me. The custom simulation script might be easier to setup import "path" for me.
Thanks for your help.
ouyang
Hi there,
My previous message might not be straight.
The point-load-cube model is shown in the attached figure in
https://groups.google.com/forum/#!topic/sfepy-devel/lBFi4LvHXZY
The Simulation started with "simple.py" works with the definition: materials = {'Load' : ({'.val' : [0.0,0.0, -100000.0]},),} equations = {'eq_1' : """dw_lin_elastic_iso.2.Omega( m.lam, m.mu, v, u ) = dw_point_load.0.Topload(Load.val, v)""",}
But Something is wrong when I follow the interactive example at: http://sfepy.org/doc-devel/tutorial.html#interactive-example-linear-elastici... with definition for the point load:
# point load is appled at vertex 4
topload= domain.create_region('Topload', 'vertex 4', 'vertex') f = Material('f', val= [[0.0],[0.0], [-100000.0]]) integral0 = Integral('i0', order=0)
t3 = Term.new('dw_point_load(f.val, v)', integral0, topload, f=f, v=v) ....
I dig into sfepy/terms/terms.py and do not see any error from the staticmethod new() in class Term for creating "dw_point_load" object.
Any suggestions would be appreciated very much.
Thanks a lot in advance.
Ouyang
On Wednesday, July 16, 2014 11:15:41 AM UTC+8, Ouyang wrote:
Hi there,
I followed the interactive example in sfepy tutorial and wrote my custom simulation script for the block domain with body force (t2 term). That works. But for the point loaded block model (t3 term), I can't make it work. The script stopped at "assert_(mat.shape[-1]==virtual.dim)" of terms_point.py with IndexError: tuple index out of range. Searching the web does not find help information for this.
The reason I write the attached custom simulation script is that I failed
doing the things below:
the "non-custom" script for a point loaded block running with "simple.py"
works and gmsh can mesh different geometries. I would like to finish those
two things together in python like "os.system('gmsh ........')", and
"os.system('simple.py examples/myscript.py") but errors to import config
in simple.py got me there since "path" settings for import in simple.py
seems a little complicated for me. The custom simulation script might be
easier to setup import "path" for me.
Thanks for your help.
ouyang
On 07/17/2014 06:21 AM, Ouyang wrote:
Hi there,
My previous message might not be straight.
The point-load-cube model is shown in the attached figure in https://groups.google.com/forum/#!topic/sfepy-devel/lBFi4LvHXZY
The Simulation started with "simple.py" works with the definition: materials = {'Load' : ({'.val' : [0.0,0.0, -100000.0]},),} equations = {'eq_1' : """dw_lin_elastic_iso.2.Omega( m.lam, m.mu, v, u ) = dw_point_load.0.Topload(Load.val, v)""",}
I have answered in the previous e-mail - just note the '.' in '.val' (special material parameter). The interactive use is under-documented in this respect.
r.
But Something is wrong when I follow the interactive example at: http://sfepy.org/doc-devel/tutorial.html#interactive-example-linear-elastici... with definition for the point load:
# point load is appled at vertex 4
topload= domain.create_region('Topload', 'vertex 4', 'vertex') f = Material('f', val= [[0.0],[0.0], [-100000.0]]) integral0 = Integral('i0', order=0)
t3 = Term.new('dw_point_load(f.val, v)', integral0, topload, f=f, v=v) ....
I dig into sfepy/terms/terms.py and do not see any error from the staticmethod new() in class Term for creating "dw_point_load" object.
Any suggestions would be appreciated very much.
Thanks a lot in advance.
Ouyang
On Wednesday, July 16, 2014 11:15:41 AM UTC+8, Ouyang wrote:
Hi there,
I followed the interactive example in sfepy tutorial and wrote my custom simulation script for the block domain with body force (t2 term). That works. But for the point loaded block model (t3 term), I can't make it work. The script stopped at "assert_(mat.shape[-1]==virtual.dim)" of terms_point.py with IndexError: tuple index out of range. Searching the web does not find help information for this.
The reason I write the attached custom simulation script is that I failed doing the things below: the "non-custom" script for a point loaded block running with "simple.py" works and gmsh can mesh different geometries. I would like to finish those two things together in python like "os.system('gmsh ........')", and "os.system('simple.py examples/myscript.py") but errors to import config in simple.py got me there since "path" settings for import in simple.py seems a little complicated for me. The custom simulation script might be easier to setup import "path" for me.
Thanks for your help.
ouyang
Hi Ouyang,
On 07/16/2014 05:15 AM, Ouyang wrote:
Hi there,
I followed the interactive example in sfepy tutorial and wrote my custom simulation script for the block domain with body force (t2 term). That works. But for the point loaded block model (t3 term), I can't make it work. The script stopped at "assert_(mat.shape[-1]==virtual.dim)" of terms_point.py with IndexError: tuple index out of range. Searching the web does not find help information for this.
It is caused by the fact, that the point load has to be given as a special material parameter (= not evaluated in quadrature points, taken "as is", see the term docstring) -> its name has to start with dot.
Try:
f = Material('f', values={'.val': [[0.0,0.0, -100000.0]]})
It works for me.
Cheers, r.
>
The reason I write the attached custom simulation script is that I failed doing the things below: the "non-custom" script for a point loaded block running with "simple.py" works and gmsh can mesh different geometries. I would like to finish those two things together in python like "os.system('gmsh ........')", and "os.system('simple.py examples/myscript.py") but errors to import config in simple.py got me there since "path" settings for import in simple.py seems a little complicated for me. The custom simulation script might be easier to setup import "path" for me.
Thanks for your help.
ouyang
Hi Robert,
hanks for your help. It works now.
I did try f = Material('f', val={'.val': [0.0,0.0, -100000.0]}) for point load and failed due to val= {...} instead of values={...}.
Based on that, I have retested the following cases both body forces and point load with three components. It seems that the argument formats of Material() are not consistency for me.
For body forces (three components), abbreviation "val=" must be used and "values=" won't work. I tested three cases: f1,f2 work and f3 fail.
f1 = Material('f', val=[[0.02], [0.01],[0.05]]) f2 = Material('f', val=[0.02, 0.01, 0.05]) f3 = Material('f', val=[[0.02, 0.01,0.05]])
For point load (three components), abbreviation "val=" won't work and "values=" must be used. As a result, f1 fail and f2,f3 work.
f1 = Material('f', values={'.val': [[0.0],[0.0], [-100000.0]]}) f2 = Material('f', values={'.val': [0.0,0.0, -100000.0]}) f3 = Material('f', values={'.val': [[0.0,0.0, -100000.0]]})
Ouyang
On Friday, July 18, 2014 6:31:14 PM UTC+8, Robert Cimrman wrote:
Hi Ouyang,
On 07/16/2014 05:15 AM, Ouyang wrote:
Hi there,
I followed the interactive example in sfepy tutorial and wrote my custom simulation script for the block domain with body force (t2 term). That works. But for the point loaded block model (t3 term), I can't make it work. The script stopped at "assert_(mat.shape[-1]==virtual.dim)" of terms_point.py with IndexError: tuple index out of range. Searching the web does not find help information for this.
It is caused by the fact, that the point load has to be given as a special material parameter (= not evaluated in quadrature points, taken "as is", see the term docstring) -> its name has to start with dot.
Try:
f = Material('f', values={'.val': [[0.0,0.0, -100000.0]]})
It works for me.
Cheers, r.
The reason I write the attached custom simulation script is that I failed doing the things below: the "non-custom" script for a point loaded block running with "simple.py" works and gmsh can mesh different geometries. I would like to finish those two things together in python like "os.system('gmsh ........')", and "os.system('simple.py examples/myscript.py") but errors to import config in simple.py got me there since "path" settings for import in simple.py seems a little complicated for me. The custom simulation script might be easier to setup import "path" for me.
Thanks for your help.
ouyang
On 07/18/2014 01:49 PM, Ouyang wrote:
Hi Robert,
hanks for your help. It works now.
I did try f = Material('f', val={'.val': [0.0,0.0, -100000.0]}) for point load and failed due to val= {...} instead of values={...}.
Here 'val' is the name you chose for the material parameter (it could be something else, like 'c' or 'const' or 'coef'...). On the other hand, values is the name of the keyword argument of Material.__init__(), that can be used to pass in material parameter names that are not valid argument names (e.g. '.val').
Based on that, I have retested the following cases both body forces and point load with three components. It seems that the argument formats of Material() are not consistency for me.
For body forces (three components), abbreviation "val=" must be used and "values=" won't work. I tested three cases: f1,f2 work and f3 fail.
f1 = Material('f', val=[[0.02], [0.01],[0.05]]) f2 = Material('f', val=[0.02, 0.01, 0.05]) f3 = Material('f', val=[[0.02, 0.01,0.05]])
For point load (three components), abbreviation "val=" won't work and "values=" must be used. As a result, f1 fail and f2,f3 work.
f1 = Material('f', values={'.val': [[0.0],[0.0], [-100000.0]]}) f2 = Material('f', values={'.val': [0.0,0.0, -100000.0]}) f3 = Material('f', values={'.val': [[0.0,0.0, -100000.0]]})
You are right that in body forces the parameter has to be a column vector, while in point loads as rows of a 2D array (each row for a point). The point load term is somewhat special.
r.
Ouyang
On Friday, July 18, 2014 6:31:14 PM UTC+8, Robert Cimrman wrote:
Hi Ouyang,
On 07/16/2014 05:15 AM, Ouyang wrote:
Hi there,
I followed the interactive example in sfepy tutorial and wrote my
custom simulation script for the block domain with body force (t2 term). That works. But for the point loaded block model (t3 term), I can't make it work. The script stopped at "assert_(mat.shape[-1]==virtual.dim)" of terms_point.py with IndexError: tuple index out of range. Searching the web does not find help information for this.
It is caused by the fact, that the point load has to be given as a special material parameter (= not evaluated in quadrature points, taken "as is", see the term docstring) -> its name has to start with dot.
Try:
f = Material('f', values={'.val': [[0.0,0.0, -100000.0]]})
It works for me.
Cheers, r.
>
The reason I write the attached custom simulation script is that I failed doing the things below: the "non-custom" script for a point loaded block running with "simple.py" works and gmsh can mesh different geometries. I would like to finish those two things together in python like "os.system('gmsh ........')", and "os.system('simple.py examples/myscript.py") but errors to import config in simple.py got me there since "path" settings for import in simple.py seems a little complicated for me. The custom simulation script might be easier to setup import "path" for me.
Thanks for your help.
ouyang
Hi Robert,
Thanks a lot for your patient explanation.
ouyang
On Friday, July 18, 2014 8:32:56 PM UTC+8, Robert Cimrman wrote:
On 07/18/2014 01:49 PM, Ouyang wrote:
Hi Robert,
hanks for your help. It works now.
I did try f = Material('f', val={'.val': [0.0,0.0, -100000.0]}) for point load and failed due to val= {...} instead of values={...}.
Here 'val' is the name you chose for the material parameter (it could be something else, like 'c' or 'const' or 'coef'...). On the other hand, values is the name of the keyword argument of Material.__init__(), that can be used to pass in material parameter names that are not valid argument names (e.g. '.val').
Based on that, I have retested the following cases both body forces and point load with three components. It seems that the argument formats of Material() are not consistency for me.
For body forces (three components), abbreviation "val=" must be used and "values=" won't work. I tested three cases: f1,f2 work and f3 fail.
f1 = Material('f', val=[[0.02], [0.01],[0.05]]) f2 = Material('f', val=[0.02, 0.01, 0.05]) f3 = Material('f', val=[[0.02, 0.01,0.05]])
For point load (three components), abbreviation "val=" won't work and "values=" must be used. As a result, f1 fail and f2,f3 work.
f1 = Material('f', values={'.val': [[0.0],[0.0], [-100000.0]]}) f2 = Material('f', values={'.val': [0.0,0.0, -100000.0]}) f3 = Material('f', values={'.val': [[0.0,0.0, -100000.0]]})
You are right that in body forces the parameter has to be a column vector, while in point loads as rows of a 2D array (each row for a point). The point load term is somewhat special.
r.
Ouyang
On Friday, July 18, 2014 6:31:14 PM UTC+8, Robert Cimrman wrote:
Hi Ouyang,
On 07/16/2014 05:15 AM, Ouyang wrote:
Hi there,
I followed the interactive example in sfepy tutorial and wrote my
custom simulation script for the block domain with body force (t2 term). That works. But for the point loaded block model (t3 term), I can't make it work. The script stopped at "assert_(mat.shape[-1]==virtual.dim)" of terms_point.py with IndexError: tuple index out of range. Searching the web does not find help information for this.
It is caused by the fact, that the point load has to be given as a special material parameter (= not evaluated in quadrature points, taken "as is", see the term docstring) -> its name has to start with dot.
Try:
f = Material('f', values={'.val': [[0.0,0.0, -100000.0]]})
It works for me.
Cheers, r.
The reason I write the attached custom simulation script is that I failed doing the things below: the "non-custom" script for a point loaded block running with "simple.py" works and gmsh can mesh different geometries. I would like to finish those two things together in python like "os.system('gmsh ........')", and "os.system('simple.py examples/myscript.py") but errors to import config in simple.py got me there since "path" settings for import in simple.py seems a little complicated for me. The custom simulation script might be easier to setup import "path" for me.
Thanks for your help.
ouyang
Hi Robert，
The attached file is basically from the interactive example of linear
elasticity and can be run by "python block_pointload_01.py". Displacements can be viewed in vtk. Can you drop several lines of code to demonstrate
how to add stress and strain output into vtk?
I did some testing between lines 85 and 105 in the attached file for that purpose but failed.
Thanks in advance.
ouyang
On Friday, July 18, 2014 8:32:56 PM UTC+8, Robert Cimrman wrote:
On 07/18/2014 01:49 PM, Ouyang wrote:
Hi Robert,
hanks for your help. It works now.
I did try f = Material('f', val={'.val': [0.0,0.0, -100000.0]}) for point load and failed due to val= {...} instead of values={...}.
Here 'val' is the name you chose for the material parameter (it could be something else, like 'c' or 'const' or 'coef'...). On the other hand, values is the name of the keyword argument of Material.__init__(), that can be used to pass in material parameter names that are not valid argument names (e.g. '.val').
Based on that, I have retested the following cases both body forces and point load with three components. It seems that the argument formats of Material() are not consistency for me.
For body forces (three components), abbreviation "val=" must be used and "values=" won't work. I tested three cases: f1,f2 work and f3 fail.
f1 = Material('f', val=[[0.02], [0.01],[0.05]]) f2 = Material('f', val=[0.02, 0.01, 0.05]) f3 = Material('f', val=[[0.02, 0.01,0.05]])
For point load (three components), abbreviation "val=" won't work and "values=" must be used. As a result, f1 fail and f2,f3 work.
f1 = Material('f', values={'.val': [[0.0],[0.0], [-100000.0]]}) f2 = Material('f', values={'.val': [0.0,0.0, -100000.0]}) f3 = Material('f', values={'.val': [[0.0,0.0, -100000.0]]})
You are right that in body forces the parameter has to be a column vector, while in point loads as rows of a 2D array (each row for a point). The point load term is somewhat special.
r.
Ouyang
On Friday, July 18, 2014 6:31:14 PM UTC+8, Robert Cimrman wrote:
Hi Ouyang,
On 07/16/2014 05:15 AM, Ouyang wrote:
Hi there,
I followed the interactive example in sfepy tutorial and wrote my
custom simulation script for the block domain with body force (t2 term). That works. But for the point loaded block model (t3 term), I can't make it work. The script stopped at "assert_(mat.shape[-1]==virtual.dim)" of terms_point.py with IndexError: tuple index out of range. Searching the web does not find help information for this.
It is caused by the fact, that the point load has to be given as a special material parameter (= not evaluated in quadrature points, taken "as is", see the term docstring) -> its name has to start with dot.
Try:
f = Material('f', values={'.val': [[0.0,0.0, -100000.0]]})
It works for me.
Cheers, r.
The reason I write the attached custom simulation script is that I failed doing the things below: the "non-custom" script for a point loaded block running with "simple.py" works and gmsh can mesh different geometries. I would like to finish those two things together in python like "os.system('gmsh ........')", and "os.system('simple.py examples/myscript.py") but errors to import config in simple.py got me there since "path" settings for import in simple.py seems a little complicated for me. The custom simulation script might be easier to setup import "path" for me.
Thanks for your help.
ouyang
On 07/28/2014 12:41 PM, Ouyang wrote:
Hi Robert，
The attached file is basically from the interactive example of linear
elasticity and can be run by "python block_pointload_01.py". Displacements can be viewed in vtk. Can you drop several lines of code to demonstrate
how to add stress and strain output into vtk?
I did some testing between lines 85 and 105 in the attached file for that purpose but failed.
Try this:
young = 2000.0 # Young's modulus [MPa]
poisson = 0.4 # Poisson's ratio
m1 = Material('m1',
values={'D': stiffness_from_youngpoisson(3, young, poisson)})
ev = pb.evaluate
strain = ev('ev_cauchy_strain.3.Omega(u)', mode='el_avg')
stress = ev('ev_cauchy_stress.3.Omega(m1.D, u)', m1=m1, mode='el_avg')
out = vec.create_output_dict()
out['cauchy_strain'] = Struct(name='output_data', mode='cell',
data=strain)
out['cauchy_stress'] = Struct(name='output_data', mode='cell',
data=stress)
pb.save_state('block_pointload_script2.vtk', out=out)
r.
Hi Robert,
It works now. Thanks a lot.
Since I learn fenics a little and run several examples, for your reference, I have a question about the "term".
Suppose we have a weak form of PDE,
\int_\Omega \Delta u \Delta v dx = \int_\Omega 6 v dx
In Fenics, this inner product, \Delta u \Delta v, can be implemented in python as
a = inner(nabla_grad(u), nabla_grad(v))*dx
it is more readable mathematically and seems any new terms can be defined here(I am not sure).
In sfepy, is there a way to do the similar thing?
Ouyang
Fenics python script for the example: ..... bc = DirichletBC(V, u0, u0_boundary)
# Define variational problem
u = TrialFunction(V) v = TestFunction(V) f = Constant(-6.0) a = inner(nabla_grad(u), nabla_grad(v))dx L = fv*dx
# Compute solution
u = Function(V) solve(a == L, u, bc) ....
On Monday, July 28, 2014 9:45:00 PM UTC+8, Robert Cimrman wrote:
On 07/28/2014 12:41 PM, Ouyang wrote:
Hi Robert，
The attached file is basically from the interactive example of
linear elasticity and can be run by "python block_pointload_01.py". Displacements can be viewed in vtk. Can you drop several lines of code to demonstrate
how to add stress and strain output into vtk?
I did some testing between lines 85 and 105 in the attached file for that purpose but failed.
Try this:
young = 2000.0 # Young's modulus [MPa]
poisson = 0.4 # Poisson's ratio
m1 = Material('m1',
values={'D': stiffness_from_youngpoisson(3, young,
poisson)})
ev = pb.evaluate
strain = ev('ev_cauchy_strain.3.Omega(u)', mode='el_avg')
stress = ev('ev_cauchy_stress.3.Omega(m1.D, u)', m1=m1, mode='el_avg')
out = vec.create_output_dict()
out['cauchy_strain'] = Struct(name='output_data', mode='cell',
data=strain)
out['cauchy_stress'] = Struct(name='output_data', mode='cell',
data=stress)
pb.save_state('block_pointload_script2.vtk', out=out)
r.
On 07/29/2014 05:38 AM, Ouyang wrote:
Hi Robert,
It works now. Thanks a lot.
Hth!
Since I learn fenics a little and run several examples, for your reference, I have a question about the "term".
Suppose we have a weak form of PDE,
\int_\Omega \Delta u \Delta v dx = \int_\Omega 6 v dx
In Fenics, this inner product, \Delta u \Delta v, can be implemented in python as
a = inner(nabla_grad(u), nabla_grad(v))*dx
it is more readable mathematically and seems any new terms can be defined here(I am not sure).
In sfepy, is there a way to do the similar thing?
No. I do not say it would be impossible to have a similar functionality, but somebody would have to think a lot and implement it.
r.
Ouyang
Fenics python script for the example: ..... bc = DirichletBC(V, u0, u0_boundary)
# Define variational problem
u = TrialFunction(V) v = TestFunction(V) f = Constant(-6.0) a = inner(nabla_grad(u), nabla_grad(v))dx L = fv*dx
# Compute solution
u = Function(V) solve(a == L, u, bc) ....
On Monday, July 28, 2014 9:45:00 PM UTC+8, Robert Cimrman wrote:
On 07/28/2014 12:41 PM, Ouyang wrote:
Hi Robert，
The attached file is basically from the interactive example of
linear elasticity and can be run by "python block_pointload_01.py". Displacements can be viewed in vtk. Can you drop several lines of code to demonstrate
how to add stress and strain output into vtk?
I did some testing between lines 85 and 105 in the attached file for that purpose but failed.
Try this:
young = 2000.0 # Young's modulus [MPa]
poisson = 0.4 # Poisson's ratio
m1 = Material('m1',
values={'D': stiffness_from_youngpoisson(3, young,
poisson)})
ev = pb.evaluate
strain = ev('ev_cauchy_strain.3.Omega(u)', mode='el_avg')
stress = ev('ev_cauchy_stress.3.Omega(m1.D, u)', m1=m1, mode='el_avg')
out = vec.create_output_dict()
out['cauchy_strain'] = Struct(name='output_data', mode='cell',
data=strain)
out['cauchy_stress'] = Struct(name='output_data', mode='cell',
data=stress)
pb.save_state('block_pointload_script2.vtk', out=out)
r.
Hi Robert,
are there any function to output just a single component of u_i (displacements) or \sigma_{ij} into vtk file?
how to output global coordinates for \sigma_{ij} in svar at qp?
Sorry. I have too many questions.
ouyang
On Tuesday, July 29, 2014 3:07:02 PM UTC+8, Robert Cimrman wrote:
On 07/29/2014 05:38 AM, Ouyang wrote:
Hi Robert,
It works now. Thanks a lot.
Hth!
Since I learn fenics a little and run several examples, for your reference, I have a question about the "term".
Suppose we have a weak form of PDE,
\int_\Omega \Delta u \Delta v dx = \int_\Omega 6 v dx
In Fenics, this inner product, \Delta u \Delta v, can be implemented in python as
a = inner(nabla_grad(u), nabla_grad(v))*dx
it is more readable mathematically and seems any new terms can be defined here(I am not sure).
In sfepy, is there a way to do the similar thing?
No. I do not say it would be impossible to have a similar functionality, but somebody would have to think a lot and implement it.
r.
Ouyang
Fenics python script for the example: ..... bc = DirichletBC(V, u0, u0_boundary)
# Define variational problem
u = TrialFunction(V) v = TestFunction(V) f = Constant(-6.0) a = inner(nabla_grad(u), nabla_grad(v))dx L = fv*dx
# Compute solution
u = Function(V) solve(a == L, u, bc) ....
On Monday, July 28, 2014 9:45:00 PM UTC+8, Robert Cimrman wrote:
On 07/28/2014 12:41 PM, Ouyang wrote:
Hi Robert，
The attached file is basically from the interactive example of
linear elasticity and can be run by "python block_pointload_01.py". Displacements can be viewed in vtk. Can you drop several lines of code to demonstrate
how to add stress and strain output into vtk?
I did some testing between lines 85 and 105 in the attached file for that purpose but failed.
Try this:
young = 2000.0 # Young's modulus [MPa]
poisson = 0.4 # Poisson's ratio
m1 = Material('m1',
values={'D': stiffness_from_youngpoisson(3, young,
poisson)})
ev = pb.evaluate
strain = ev('ev_cauchy_strain.3.Omega(u)', mode='el_avg')
stress = ev('ev_cauchy_stress.3.Omega(m1.D, u)', m1=m1,
mode='el_avg')
out = vec.create_output_dict()
out['cauchy_strain'] = Struct(name='output_data', mode='cell',
data=strain)
out['cauchy_stress'] = Struct(name='output_data', mode='cell',
data=stress)
pb.save_state('block_pointload_script2.vtk', out=out)
r.
On 07/31/2014 12:12 PM, Ouyang wrote:
Hi Robert,
I guess you do not need to interpolate the stresses into a field variable - you can work with the element-averaged values, no?
BTW. a function for computing principal stresses would be a nice addition to sfepy.mechanics.tensors.
Not per se, but it's really easy - just add that to the output dictionary:
from sfepy.mechanics.tensors import get_full_indices
s2f = nm.array(get_full_indices(3))
ii = s2f[0, 1]
out['cauchy_stress_01'] = Struct(name='output_data', mode='cell',
data=stress[..., ii:ii+1, :])
out['u_0'] = Struct(name='output_data', mode='vertex',
data = out['u'].data[:, 0:1])
You mean how to get the physical coordinates of all the qp in the whole domain? Or something else?
Sorry. I have too many questions.
No problem :) But you will have to explain more what are you trying to do.
r.
On Thursday, July 31, 2014 6:56:46 PM UTC+8, Robert Cimrman wrote: >
On 07/31/2014 12:12 PM, Ouyang wrote:
Hi Robert,
I guess you do not need to interpolate the stresses into a field variable
you can work with the element-averaged values, no?
I want to evaluate the stresses at center of the block which point load is applied. For brittle materials, rock, we concern its tensile stress (normally minimum stress defined in rock mechanics) at the block center to see if it reaches its tensile strength. So, principal stresses is normally easy to be used in a failure criteria in engineering. I guess the current element types available in sfepy are constant-strain or constant-stress in a element. So my understanding is that the constant-strain or constant-stress is element-averaged values.
BTW. a function for computing principal stresses would be a nice addition to sfepy.mechanics.tensors.
I am not sure I am good enough for this addition. I can try it when I have time.
>
Not per se, but it's really easy - just add that to the output dictionary:
from sfepy.mechanics.tensors import get_full_indices
s2f = nm.array(get_full_indices(3))
ii = s2f[0, 1]
out['cauchy_stress_01'] = Struct(name='output_data', mode='cell',
data=stress[..., ii:ii+1, :])
out['u_0'] = Struct(name='output_data', mode='vertex',
data = out['u'].data[:, 0:1])
Thanks for this. I will test it.
You mean how to get the physical coordinates of all the qp in the whole domain? Or something else?
Yes, the physical coordinates associated with the stresses/strains/displacements in the whole domain. It could be the qp or some other points. If so, probably more things we can do in the postprocessing.
Sorry. I have too many questions.
No problem :) But you will have to explain more what are you trying to do.
r.
Thanks again. ouyang
On 07/31/2014 02:27 PM, Ouyang wrote: > >
On Thursday, July 31, 2014 6:56:46 PM UTC+8, Robert Cimrman wrote: >
On 07/31/2014 12:12 PM, Ouyang wrote:
Hi Robert,
1. how to output principal stresses in sfepy? I can calculate
principal stresses based on \sigma_{ij} in svar below. But I guess it is not straightforward for the task in my question 2. ... sfield = Field.from_args('stress', nm.float64, (6,), pb.domain.regions['Omega']) svar = FieldVariable('sigma', 'parameter', sfield,primary_var_name='(set-to-None)') svar.set_data_from_qp(stress, ivn) ...
I guess you do not need to interpolate the stresses into a field variable
you can work with the element-averaged values, no?
I want to evaluate the stresses at center of the block which point load is applied. For brittle materials, rock, we concern its tensile stress (normally minimum stress defined in rock mechanics) at the block center to see if it reaches its tensile strength. So, principal stresses is normally easy to be used in a failure criteria in engineering.
I see. So set_data_from_qp() might be good for you. What it actually does is that it integrates a quantity (given in qp) in all elements around a vertex and divides it by the volume of those elements. So it's a volume-weighted average.
I guess the current element types available in sfepy are constant-strain or constant-stress in a element. So my understanding is that the constant-strain or constant-stress is element-averaged values.
Not only constant-strain - this depends on the field approximation order. If you set the displacement field order to something > 1, even for tetrahedra the train in an element will not be constant.
>
BTW. a function for computing principal stresses would be a nice addition to sfepy.mechanics.tensors.
I am not sure I am good enough for this addition. I can try it when I have time.
Thanks! (Are you interested just in principal stress values, or also in the directions?)
>
2. are there any function to output just a single component of
u_i
(displacements) or \sigma_{ij} into vtk file?
Not per se, but it's really easy - just add that to the output dictionary:
from sfepy.mechanics.tensors import get_full_indices
s2f = nm.array(get_full_indices(3))
ii = s2f[0, 1]
out['cauchy_stress_01'] = Struct(name='output_data', mode='cell',
data=stress[..., ii:ii+1, :])
out['u_0'] = Struct(name='output_data', mode='vertex',
data = out['u'].data[:, 0:1])
Thanks for this. I will test it.
3. how to output global coordinates for
\sigma_{ij} in svar at qp?
You mean how to get the physical coordinates of all the qp in the whole domain? Or something else?
Yes, the physical coordinates associated with the stresses/strains/displacements in the whole domain. It could be the qp or some other points. If so, probably more things we can do in the postprocessing.
Try Term.get_physical_qps() with your term instances, or something like:
from sfepy.discrete.common.mappings import get_physical_qps qps = get_physical_qps(omega, integral) coors = qps.get_merged_values()
The QP coordinates are given in the order of the mesh elements.
For the coordinates of the field nodes (for the lagrange basis), you can use FEField.get_coor() with your field instances (e.g. u.field.get_coor()). For linear fields covering the whole domain (as in your case), the coordinates are the same as the mesh vertices (domain.mesh.coors or mesh.coors).
Does that help?
r.
On Thursday, July 31, 2014 8:49:23 PM UTC+8, Robert Cimrman wrote: >
On 07/31/2014 02:27 PM, Ouyang wrote:
On Thursday, July 31, 2014 6:56:46 PM UTC+8, Robert Cimrman wrote:
On 07/31/2014 12:12 PM, Ouyang wrote:
Hi Robert,
1. how to output principal stresses in sfepy? I can calculate
principal stresses based on \sigma_{ij} in svar below. But I guess it is not straightforward for the task in my question 2. ... sfield = Field.from_args('stress', nm.float64, (6,), pb.domain.regions['Omega']) svar = FieldVariable('sigma', 'parameter', sfield,primary_var_name='(set-to-None)') svar.set_data_from_qp(stress, ivn) ...
I guess you do not need to interpolate the stresses into a field variable
you can work with the element-averaged values, no?
I want to evaluate the stresses at center of the block which point load is applied. For brittle materials, rock, we concern its tensile stress (normally minimum stress defined in rock mechanics) at the block center to see if it reaches its tensile strength. So, principal stresses is normally easy to be used in a failure criteria in engineering.
I see. So set_data_from_qp() might be good for you. What it actually does is that it integrates a quantity (given in qp) in all elements around a vertex and divides it by the volume of those elements. So it's a volume-weighted average.
I guess the current element types available in sfepy are constant-strain or constant-stress in a element. So my understanding is that the constant-strain or constant-stress is element-averaged values.
Not only constant-strain - this depends on the field approximation order. If you set the displacement field order to something > 1, even for tetrahedra the train in an element will not be constant.
BTW. a function for computing principal stresses would be a nice addition to sfepy.mechanics.tensors.
I am not sure I am good enough for this addition. I can try it when I have time.
Thanks! (Are you interested just in principal stress values, or also in the directions?)
Yes. both magnitudes and directions. I used to use "eigenvalues and eigenvectors" http://docs.scipy.org/doc/scipy/reference/tutorial/linalg.html#eigenvalues-a... from scipy:
import numpy as np from scipy import linalg A = np.array([[1,5,2],[2,4,1],[3,6,2]]) la,v = linalg.eig(A)
la are the magnitudes of principal stresses and v holds its directions.
Does this "fit" somewhere in tensors.py of sfepy?
I guess not that easy.
2. are
there any function to output just a single component of u_i
(displacements) or \sigma_{ij} into vtk file?
Not per se, but it's really easy - just add that to the output dictionary:
from sfepy.mechanics.tensors import get_full_indices
s2f = nm.array(get_full_indices(3))
ii = s2f[0, 1]
out['cauchy_stress_01'] = Struct(name='output_data', mode='cell',
data=stress[..., ii:ii+1, :])
out['u_0'] = Struct(name='output_data', mode='vertex',
data = out['u'].data[:, 0:1])
Thanks for this. I will test it.
3. how to output global coordinates for
\sigma_{ij} in svar at qp?
You mean how to get the physical coordinates of all the qp in the whole domain? Or something else?
Yes, the physical coordinates associated with the stresses/strains/displacements in the whole domain. It could be the qp or some other points. If so, probably more things we can do in the postprocessing.
Try Term.get_physical_qps() with your term instances, or something like:
from sfepy.discrete.common.mappings import get_physical_qps qps = get_physical_qps(omega, integral) coors = qps.get_merged_values()
The QP coordinates are given in the order of the mesh elements.
For the coordinates of the field nodes (for the lagrange basis), you can use FEField.get_coor() with your field instances (e.g. u.field.get_coor()). For linear fields covering the whole domain (as in your case), the coordinates are the same as the mesh vertices (domain.mesh.coors or mesh.coors).
Does that help?
r.
Yes. that helps a lot. I will test it.
Thanks.
ouyang
On 07/31/2014 05:53 PM, Ouyang wrote:
BTW. a function for computing principal stresses would be a nice addition to sfepy.mechanics.tensors.
I am not sure I am good enough for this addition. I can try it when I have time.
Thanks! (Are you interested just in principal stress values, or also in the directions?)
Yes. both magnitudes and directions. I used to use "eigenvalues and eigenvectors" http://docs.scipy.org/doc/scipy/reference/tutorial/linalg.html#eigenvalues-a... from scipy:
import numpy as np from scipy import linalg A = np.array([[1,5,2],[2,4,1],[3,6,2]]) la,v = linalg.eig(A)
la are the magnitudes of principal stresses and v holds its directions.
Does this "fit" somewhere in tensors.py of sfepy?
I guess not that easy.
I was thinking in terms of a dedicated function that would compute the eigenvalues and eigenvectors of a symmetric 3x3 or 2x2 matrix in a vectorized way, so that an array of stress tensors could be passed in. That would be handy, but it is not a _must have now_ feature :) It might fit better somewhere into sfepy.linalg.
Try Term.get_physical_qps() with your term instances, or something like:
from sfepy.discrete.common.mappings import get_physical_qps qps = get_physical_qps(omega, integral) coors = qps.get_merged_values()
The QP coordinates are given in the order of the mesh elements.
For the coordinates of the field nodes (for the lagrange basis), you can use FEField.get_coor() with your field instances (e.g. u.field.get_coor()). For linear fields covering the whole domain (as in your case), the coordinates are the same as the mesh vertices (domain.mesh.coors or mesh.coors).
Does that help?
r.
Yes. that helps a lot. I will test it.
OK, let us know how it goes.
r.
On Thursday, July 31, 2014 11:59:06 PM UTC+8, Robert Cimrman wrote: >
On 07/31/2014 05:53 PM, Ouyang wrote:
BTW. a function for computing principal stresses would be a nice addition to sfepy.mechanics.tensors.
I am not sure I am good enough for this addition. I can try it when I have time.
Thanks! (Are you interested just in principal stress values, or also in the directions?)
Yes. both magnitudes and directions. I used to use "eigenvalues and eigenvectors" http://docs.scipy.org/doc/scipy/reference/tutorial/linalg.html#eigenvalues-... from scipy:
import numpy as np from scipy import linalg A = np.array([[1,5,2],[2,4,1],[3,6,2]]) la,v = linalg.eig(A)
la are the magnitudes of principal stresses and v holds its directions.
Does this "fit" somewhere in tensors.py of sfepy?
I guess not that easy.
I was thinking in terms of a dedicated function that would compute the eigenvalues and eigenvectors of a symmetric 3x3 or 2x2 matrix in a vectorized way, so that an array of stress tensors could be passed in. That would be handy, but it is not a _must have now_ feature :) It might fit better somewhere into sfepy.linalg.
I am still in the early stage of learning sfepy and too many new things for me :)).
Try Term.get_physical_qps() with your term instances, or something like:
from sfepy.discrete.common.mappings import get_physical_qps qps = get_physical_qps(omega, integral) coors = qps.get_merged_values()
The QP coordinates are given in the order of the mesh elements.
For the coordinates of the field nodes (for the lagrange basis), you can use FEField.get_coor() with your field instances (e.g. u.field.get_coor()). For linear fields covering the whole domain (as in your case), the coordinates are the same as the mesh vertices (domain.mesh.coors or mesh.coors).
Does that help?
r.
Yes. that helps a lot. I will test it.
OK, let us know how it goes.
r.
I have tested the scripts in your message. All works but only one thing: in the figure below, most left cube in the figure below is displacement u_3 only and most right cube is stress \sigma_33. u_3's contour looks good but \sigma_33 's contour jumps from element to element. Do you know why?
Ouyang
https://lh5.googleusercontent.com/-LpsBQHDxanQ/U9tx3S6sDxI/AAAAAAAAC4w/kZCte...
On 08/01/2014 01:02 PM, Ouyang wrote: > >
On Thursday, July 31, 2014 11:59:06 PM UTC+8, Robert Cimrman wrote: >
On 07/31/2014 05:53 PM, Ouyang wrote:
BTW. a function for computing principal stresses would be a nice addition to sfepy.mechanics.tensors.
I am not sure I am good enough for this addition. I can try it when I have time.
Thanks! (Are you interested just in principal stress values, or also in the directions?)
Yes. both magnitudes and directions. I used to use "eigenvalues and eigenvectors" http://docs.scipy.org/doc/scipy/reference/tutorial/linalg.html#eigenvalues-... from scipy:
import numpy as np from scipy import linalg A = np.array([[1,5,2],[2,4,1],[3,6,2]]) la,v = linalg.eig(A)
la are the magnitudes of principal stresses and v holds its directions.
Does this "fit" somewhere in tensors.py of sfepy?
I guess not that easy.
I was thinking in terms of a dedicated function that would compute the eigenvalues and eigenvectors of a symmetric 3x3 or 2x2 matrix in a vectorized way, so that an array of stress tensors could be passed in. That would be handy, but it is not a _must have now_ feature :) It might fit better somewhere into sfepy.linalg.
I am still in the early stage of learning sfepy and too many new things for me :)).
Yes, the basic use is quite simple, but after you get into gory details... :)
Try Term.get_physical_qps() with your term instances, or something like: >
from sfepy.discrete.common.mappings import get_physical_qps qps = get_physical_qps(omega, integral) coors = qps.get_merged_values()
The QP coordinates are given in the order of the mesh elements.
For the coordinates of the field nodes (for the lagrange basis), you can use FEField.get_coor() with your field instances (e.g. u.field.get_coor()). For linear fields covering the whole domain (as in your case), the coordinates are the same as the mesh vertices (domain.mesh.coors or mesh.coors).
Does that help?
r.
Yes. that helps a lot. I will test it.
OK, let us know how it goes.
r.
I have tested the scripts in your message. All works but only one thing: in the figure below, most left cube in the figure below is displacement u_3 only and most right cube is stress \sigma_33. u_3's contour looks good but \sigma_33 's contour jumps from element to element. Do you know why?
Ouyang
https://lh5.googleusercontent.com/-LpsBQHDxanQ/U9tx3S6sDxI/AAAAAAAAC4w/kZCte...
This is a well-known problem for constant-stress elements, see for example [1], where some solution is proposed (just a quick google search returned that, I did not read it).
Try refining the mesh around that corner and/or use some averaging.
r.
[1] http://web.mit.edu/kjb/www/Principal_Publications/A_stress_improvement_proce...
I forgot to attach the file for your reference in the previous message. ouyang
On Tuesday, July 29, 2014 3:07:02 PM UTC+8, Robert Cimrman wrote:
On 07/29/2014 05:38 AM, Ouyang wrote:
Hi Robert,
It works now. Thanks a lot.
Hth!
Since I learn fenics a little and run several examples, for your reference, I have a question about the "term".
Suppose we have a weak form of PDE,
\int_\Omega \Delta u \Delta v dx = \int_\Omega 6 v dx
In Fenics, this inner product, \Delta u \Delta v, can be implemented in python as
a = inner(nabla_grad(u), nabla_grad(v))*dx
it is more readable mathematically and seems any new terms can be defined here(I am not sure).
In sfepy, is there a way to do the similar thing?
No. I do not say it would be impossible to have a similar functionality, but somebody would have to think a lot and implement it.
r.
Ouyang
Fenics python script for the example: ..... bc = DirichletBC(V, u0, u0_boundary)
# Define variational problem
u = TrialFunction(V) v = TestFunction(V) f = Constant(-6.0) a = inner(nabla_grad(u), nabla_grad(v))dx L = fv*dx
# Compute solution
u = Function(V) solve(a == L, u, bc) ....
On Monday, July 28, 2014 9:45:00 PM UTC+8, Robert Cimrman wrote:
On 07/28/2014 12:41 PM, Ouyang wrote:
Hi Robert，
The attached file is basically from the interactive example of
linear elasticity and can be run by "python block_pointload_01.py". Displacements can be viewed in vtk. Can you drop several lines of code to demonstrate
how to add stress and strain output into vtk?
I did some testing between lines 85 and 105 in the attached file for that purpose but failed.
Try this:
young = 2000.0 # Young's modulus [MPa]
poisson = 0.4 # Poisson's ratio
m1 = Material('m1',
values={'D': stiffness_from_youngpoisson(3, young,
poisson)})
ev = pb.evaluate
strain = ev('ev_cauchy_strain.3.Omega(u)', mode='el_avg')
stress = ev('ev_cauchy_stress.3.Omega(m1.D, u)', m1=m1,
mode='el_avg')
out = vec.create_output_dict()
out['cauchy_strain'] = Struct(name='output_data', mode='cell',
data=strain)
out['cauchy_stress'] = Struct(name='output_data', mode='cell',
data=stress)
pb.save_state('block_pointload_script2.vtk', out=out)
r.