Skip to content

Commit

Permalink
Merge branch 'reorganize-examples'
Browse files Browse the repository at this point in the history
  • Loading branch information
@rc
rc committed Mar 27, 2015
2 parents ebee3c3 + fefa790 commit fe708a4
Show file tree
Hide file tree
Showing 34 changed files with 168 additions and 250 deletions.
4 changes: 2 additions & 2 deletions doc/tutorial.rst
View file
Expand Up @@ -571,8 +571,8 @@ Complete Example as a Script
^^^^^^^^^^^^^^^^^^^^^^^^^^^^

The source code: :download:`linear_elasticity.py
<../examples/standalone/interactive/linear_elasticity.py>`. It should be
<../examples/linear_elasticity/linear_elastic_interactive.py>`. It should be
run from the *SfePy* source directory so that it finds the mesh file.

.. literalinclude:: ../examples/standalone/interactive/linear_elasticity.py
.. literalinclude:: ../examples/linear_elasticity/linear_elastic_interactive.py
:linenos:
0 ...e/interactive/laplace_shifted_periodic.py → ...les/diffusion/laplace_shifted_periodic.py
View file
File renamed without changes.
205 changes: 84 additions & 121 deletions ...e/homogenized_elasticity/rs_correctors.py → examples/homogenization/rs_correctors.py 100644 → 100755
View file
@@ -1,41 +1,45 @@
# c: 05.05.2008, r: 05.05.2008
#!/usr/bin/env python
"""
Compute homogenized elastic coefficients for a given microstructure.
"""
from optparse import OptionParser
import sys
sys.path.append( '.' )
sys.path.append('.')

import numpy as nm

from sfepy import data_dir
import sfepy.discrete.fem.periodic as per
from sfepy.homogenization.utils import define_box_regions

# c: 05.05.2008, r: 05.05.2008
def define_regions( filename ):
"""Define various subdomain for a given mesh file. This function is called
below."""
def define_regions(filename):
"""
Define various subdomains for a given mesh file.
"""
regions = {}
dim = 2

regions['Y'] = 'all'

eog = 'cells of group %d'
if filename.find( 'osteonT1' ) >= 0:
if filename.find('osteonT1') >= 0:
mat_ids = [11, 39, 6, 8, 27, 28, 9, 2, 4, 14, 12, 17, 45, 28, 15]
regions['Ym'] = ' +c '.join( (eog % im) for im in mat_ids )
regions['Ym'] = ' +c '.join((eog % im) for im in mat_ids)
wx = 0.865
wy = 0.499

regions['Yc'] = 'r.Y -c r.Ym'

# Sides and corners.
regions.update( define_box_regions( 2, (wx, wy) ) )
regions.update(define_box_regions(2, (wx, wy)))

return dim, regions, mat_ids

def get_pars(ts, coor, mode=None, term=None, group_indx=None,
mat_ids = [], **kwargs):
"""Define material parameters:
$D_ijkl$ (elasticity),
in a given region."""
"""
Define material parameters: :math:`D_ijkl` (elasticity), in a given region.
"""
if mode == 'qp':
dim = coor.shape[1]
sym = (dim + 1) * dim / 2
Expand All @@ -48,9 +52,9 @@ def get_pars(ts, coor, mode=None, term=None, group_indx=None,
# in 1e+10 [Pa]
lam = 1.7
mu = 0.3
o = nm.array( [1.] * dim + [0.] * (sym - dim), dtype = nm.float64 )
oot = nm.outer( o, o )
out['D'] = lam * oot + mu * nm.diag( o + 1.0 )
o = nm.array([1.] * dim + [0.] * (sym - dim), dtype = nm.float64)
oot = nm.outer(o, o)
out['D'] = lam * oot + mu * nm.diag(o + 1.0)

for key, val in out.iteritems():
out[key] = nm.tile(val, (coor.shape[0], 1, 1))
Expand All @@ -63,12 +67,12 @@ def get_pars(ts, coor, mode=None, term=None, group_indx=None,

##
# Mesh file.
filename_mesh = 'meshes/osteonT1_11.mesh'
filename_mesh = data_dir + '/meshes/2d/special/osteonT1_11.mesh'

##
# Define regions (subdomains, boundaries) - $Y$, $Y_i$, ...
# depending on a mesh used.
dim, regions, mat_ids = define_regions( filename_mesh )
dim, regions, mat_ids = define_regions(filename_mesh)

functions = {
'get_pars' : (lambda ts, coors, **kwargs:
Expand All @@ -83,12 +87,8 @@ def get_pars(ts, coor, mode=None, term=None, group_indx=None,
##
# Define fields: 'displacement' in $Y$,
# 'pressure_m' in $Y_m$.
field_1 = {
'name' : 'displacement',
'dtype' : nm.float64,
'shape' : dim,
'region' : 'Y',
'approx_order' : 1,
fields = {
'displacement' : ('real', dim, 'Y', 1),
}

##
Expand All @@ -105,38 +105,22 @@ def get_pars(ts, coor, mode=None, term=None, group_indx=None,
##
# Periodic boundary conditions.
if dim == 3:
epbc_10 = {
'name' : 'periodic_x',
'region' : ['Left', 'Right'],
'dofs' : {'uc.all' : 'uc.all'},
'match' : 'match_x_plane',
}
epbc_11 = {
'name' : 'periodic_y',
'region' : ['Near', 'Far'],
'dofs' : {'uc.all' : 'uc.all'},
'match' : 'match_y_plane',
}
epbc_12 = {
'name' : 'periodic_z',
'region' : ['Top', 'Bottom'],
'dofs' : {'uc.all' : 'uc.all'},
'match' : 'match_z_plane',
epbcs = {
'periodic_x' : (['Left', 'Right'], {'uc.all' : 'uc.all'},
'match_x_plane'),
'periodic_y' : (['Near', 'Far'], {'uc.all' : 'uc.all'},
'match_y_plane'),
'periodic_z' : (['Top', 'Bottom'], {'uc.all' : 'uc.all'},
'match_z_plane'),
}
else:
epbc_10 = {
'name' : 'periodic_x',
'region' : ['Left', 'Right'],
'dofs' : {'uc.all' : 'uc.all'},
'match' : 'match_y_line',
epbcs = {
'periodic_x' : (['Left', 'Right'], {'uc.all' : 'uc.all'},
'match_y_line'),
'periodic_y' : (['Bottom', 'Top'], {'uc.all' : 'uc.all'},
'match_x_line'),
}
epbc_11 = {
'name' : 'periodic_y',
'region' : ['Top', 'Bottom'],
'dofs' : {'uc.all' : 'uc.all'},
'match' : 'match_x_line',
}


##
# Dirichlet boundary conditions.
ebcs = {
Expand All @@ -145,16 +129,14 @@ def get_pars(ts, coor, mode=None, term=None, group_indx=None,

##
# Material defining constitutive parameters of the microproblem.
material_1 = {
'name' : 'm',
'function' : 'get_pars',
materials = {
'm' : 'get_pars',
}

##
# Numerical quadratures for volume (i3 - order 3) integral terms.
integral_1 = {
'name' : 'i3',
'order' : 3,
integrals = {
'i3' : 3,
}

##
Expand All @@ -169,7 +151,7 @@ def set_elastic(variables, ir, ic, mode, pis, corrs_rs):
coefs = {
'E' : {
'requires' : ['pis', 'corrs_rs'],
'expression' : 'dw_lin_elastic.i3.Y( m.D, Pi1, Pi2 )',
'expression' : 'dw_lin_elastic.i3.Y(m.D, Pi1, Pi2)',
'set_variables' : set_elastic,
},
}
Expand All @@ -186,8 +168,8 @@ def set_elastic(variables, ir, ic, mode, pis, corrs_rs):
'save_variables' : ['uc'],
'ebcs' : ['fixed_u'],
'epbcs' : all_periodic,
'equations' : {'eq' : """dw_lin_elastic.i3.Y( m.D, vc, uc )
= - dw_lin_elastic.i3.Y( m.D, vc, Pi )"""},
'equations' : {'eq' : """dw_lin_elastic.i3.Y(m.D, vc, uc)
= - dw_lin_elastic.i3.Y(m.D, vc, Pi)"""},
'set_variables' : [('Pi', 'pis', 'uc')],
'save_name' : 'corrs_elastic',
'is_linear' : True,
Expand All @@ -196,38 +178,22 @@ def set_elastic(variables, ir, ic, mode, pis, corrs_rs):

##
# Solvers.
solver_0 = {
'name' : 'ls',
'kind' : 'ls.scipy_direct', # Direct solver.
}

solver_1 = {
'name' : 'newton',
'kind' : 'nls.newton',

'i_max' : 2,
'eps_a' : 1e-8,
'eps_r' : 1e-2,
'macheps' : 1e-16,
'lin_red' : 1e-2, # Linear system error < (eps_a * lin_red).
'ls_red' : 0.1,
'ls_red_warp' : 0.001,
'ls_on' : 0.99999,
'ls_min' : 1e-5,
'check' : 0,
'delta' : 1e-6,
'problem' : 'nonlinear', # 'nonlinear' or 'linear' (ignore i_max)
solvers = {
'ls' : ('ls.scipy_direct', {}),
'newton' : ('nls.newton', {
'i_max' : 1,
'eps_a' : 1e-8,
'eps_r' : 1e-2,
})
}

############################################
# Mini-application below, computing the homogenized elastic coefficients.
usage = """%prog [options]"""
usage = '%prog [options]\n' + __doc__.rstrip()
help = {
'no_pauses' : 'do not make pauses',
}

##
# c: 05.05.2008, r: 28.11.2008
def main():
import os
from sfepy.base.base import spause, output
Expand All @@ -245,93 +211,90 @@ def main():
def spause(*args):
output(*args)

nm.set_printoptions( precision = 3 )
nm.set_printoptions(precision=3)

spause( r""">>>
spause(r""">>>
First, this file will be read in place of an input
(problem description) file.
Press 'q' to quit the example, press any other key to continue...""" )
Press 'q' to quit the example, press any other key to continue...""")
required, other = get_standard_keywords()
required.remove( 'equations' )
required.remove('equations')
# Use this file as the input file.
conf = ProblemConf.from_file( __file__, required, other )
conf = ProblemConf.from_file(__file__, required, other)
print conf.to_dict().keys()
spause( r""">>>
spause(r""">>>
...the read input as a dict (keys only for brevity).
['q'/other key to quit/continue...]""" )
['q'/other key to quit/continue...]""")

spause( r""">>>
spause(r""">>>
Now the input will be used to create a Problem instance.
['q'/other key to quit/continue...]""" )
['q'/other key to quit/continue...]""")
problem = Problem.from_conf(conf, init_equations=False)
# The homogenization mini-apps need the output_dir.
output_dir = os.path.join(os.path.split(__file__)[0], 'output')
if not os.path.exists(output_dir):
os.makedirs(output_dir)
output_dir = ''
problem.output_dir = output_dir
print problem
spause( r""">>>
spause(r""">>>
...the Problem instance.
['q'/other key to quit/continue...]""" )
['q'/other key to quit/continue...]""")


spause( r""">>>
spause(r""">>>
The homogenized elastic coefficient $E_{ijkl}$ is expressed
using $\Pi$ operators, computed now. In fact, those operators are permuted
coordinates of the mesh nodes.
['q'/other key to quit/continue...]""" )
['q'/other key to quit/continue...]""")
req = conf.requirements['pis']
mini_app = cb.ShapeDimDim( 'pis', problem, req )
mini_app = cb.ShapeDimDim('pis', problem, req)
pis = mini_app()
print pis
spause( r""">>>
spause(r""">>>
...the $\Pi$ operators.
['q'/other key to quit/continue...]""" )
['q'/other key to quit/continue...]""")

spause( r""">>>
spause(r""">>>
Next, $E_{ijkl}$ needs so called steady state correctors $\bar{\omega}^{rs}$,
computed now.
['q'/other key to quit/continue...]""" )
['q'/other key to quit/continue...]""")
req = conf.requirements['corrs_rs']

save_name = req.get( 'save_name', '' )
name = os.path.join( output_dir, save_name )
save_name = req.get('save_name', '')
name = os.path.join(output_dir, save_name)

mini_app = cb.CorrDimDim('steady rs correctors', problem, req)
mini_app.setup_output(save_format='vtk',
file_per_var=False)
corrs_rs = mini_app( data = {'pis': pis} )
corrs_rs = mini_app(data={'pis': pis})
print corrs_rs
spause( r""">>>
spause(r""">>>
...the $\bar{\omega}^{rs}$ correctors.
The results are saved in: %s.%s
Try to display them with:
python postproc.py %s.%s
['q'/other key to quit/continue...]""" % (2 * (name, problem.output_format)) )
['q'/other key to quit/continue...]""" % (2 * (name, problem.output_format)))

spause( r""">>>
spause(r""">>>
Then the volume of the domain is needed.
['q'/other key to quit/continue...]""" )
volume = problem.evaluate('d_volume.i3.Y( uc )')
['q'/other key to quit/continue...]""")
volume = problem.evaluate('d_volume.i3.Y(uc)')
print volume

spause( r""">>>
spause(r""">>>
...the volume.
['q'/other key to quit/continue...]""" )
['q'/other key to quit/continue...]""")

spause( r""">>>
spause(r""">>>
Finally, $E_{ijkl}$ can be computed.
['q'/other key to quit/continue...]""" )
['q'/other key to quit/continue...]""")
mini_app = cb.CoefSymSym('homogenized elastic tensor',
problem, conf.coefs['E'])
c_e = mini_app(volume, data={'pis': pis, 'corrs_rs' : corrs_rs})
print r""">>>
The homogenized elastic coefficient $E_{ijkl}$, symmetric storage
with rows, columns in 11, 22, 12 ordering:"""
print c_e

if __name__ == '__main__':
main()

0 comments on commit fe708a4

Please sign in to comment.