On 04/21/10 00:43, Dominique wrote:
On Apr 20, 2:43 pm, Robert
Cimrman<cimr...(a)ntc.zcu.cz> wrote:
> Hi Dominique,
> On 04/20/10 15:29, Dominique wrote:
>> Hi Robert,
>>> From a given configuration file, I
think I am able to extract the
>> Jacobian (called the "tangent matrix" in the code) and the right-hand
>> side as follows:
>> def get_problem(fName=__file__):
>> """
>> Read the configuration specified in the file specified in argument
>> (the present file by default) and return an initialized problem
>> constructed from it along with an initial state vector.
>> """
>> from sfepy.base.conf import ProblemConf, get_standard_keywords
>> from sfepy.fem import eval_term_op, ProblemDefinition
>> required, other =
get_standard_keywords()
>> # Use this file as the input file
to gather configuration.
>> conf = ProblemConf.from_file(fName, required, other)
>> # Create a ProblemDefinition
instance from current configuration.
>> problem = ProblemDefinition.from_conf(conf,
>> init_variables=True,
>> init_equations=True)
> There is a shortcut for this:
> problem =
ProblemDefinition.from_conf_file(fName,required=required, other=other)
>> state =
problem.create_state_vector()
>> problem.init_variables(state)
>> problem.time_update(None)
>> return (problem, state)
>> def get_jacobian(problem, state=None):
>> """
>> Return the Jacobian matrix of the system in SciPy CSR format.
>> Input arguments are as returned by `get_problem()`.
>> """
>> # Instantiate a problem evaluator.
>> evaluator = problem.get_evaluator(mtx=problem.mtx_a)
>> if state is None:
>> state = problem.create_state_vector()
>> jac = evaluator.eval_tangent_matrix(state)
>> return jac
>> Am I missing anything here (e.g., I'm
not sure the BCs and ICs are
>> taken into account)?
> This should be cared for by the
problem.time_update() call.
>> Now I'd like to extract the
contribution of a particular equation to
>> this Jacobian. How can I do that? This is what I've tried but it seems
>> to return the whole Jacobian:
>> def get_equation_jacobian(problem,
equation, state=None):
>> """
>> Return the Jacobian matrix of a specific equation in SciPy CSR
>> format.
>> """
>> from sfepy.fem.evaluate import assemble_matrix
>> evaluator = problem.get_evaluator(mtx=problem.mtx_a)
>> if state is None:
>> state = problem.create_state_vector()
>> problem.variables.data_from_state(state)
>> jac = problem.mtx_a
>> jac = assemble_matrix(jac, equation,
>> problem.variables, problem.materials,
>> **evaluator.data)
>> return jac
>> In fact, this function returns None, but
problem.mtx_a contains the
>> entire Jacobian. I suppose I should specify a different matrix when I
>> instantiate the evaluator, but where can I discover its size?
> The sizes of the various components can be
obtained by:
> (Pdb) print problem.variables.di
> DofInfo:dof_info
> details:
> {'p': Struct:field_var_dof_details, 'u':
Struct:field_var_dof_details}
> indx:
> {'p': slice(2484, 3312, None), 'u': slice(0, 2484, None)}
> n_dof:
> {'p': 828, 'u': 2484}
> name:
> dof_info
> names:
> ['u', 'p']
> ptr:
> [ 0 2484 3312]
> problem.variables.di is the full DOF info
(including BC)
> problem.variables.adi is the stripped DOF info (without the fixed DOFs)
> Using the 'indx', you can index a
submatrix.
> Look also at assemble_by_blocks() in
sfepy/fem/evaluate.py.
> Does that help?
Yes, although it seems I'm having more
success with eval_term_op as
follows:
def get_jacobian_blocks(problem, state=None):
"""
Return the blocks of the Jacobian matrix as a dictionary
of sparse matrices in SciPy CSR format.
Input arguments are as returned by `get_problem()`.
"""
from sfepy.fem.evaluate import eval_term_op
if state is None:
state = problem.create_state_vector()
blocks = {}
for key in problem.conf.equations.keys():
equation = problem.conf.equations[key]
blocks[key] = eval_term_op(state, equation, problem,
dw_mode='matrix',
tangent_matrix=None)
return blocks
The difficulty is that all blocks have the size
of the full Jacobian.
I guess I could extract the actual blocks that I need using the
dimensions given by problem.variables.di but I suspect there's a
better way...
Yes, your way is ok. It could be made easier though - there could be a helper
function directly in sfepy.
Maybe an example of input for assemble_by_blocks() would help too:
variables = {
'u' : ('unknown field', 'displacement_Y2', 0),
'v' : ('test field', 'displacement_Y2', 'u'),
'phi' : ('unknown field', 'potential_Y2', 0),
'psi' : ('test field', 'potential_Y2', 'phi'),
}
equations = {
'K,v,u' :
"""dw_lin_elastic.i1.Y2( inclusion.C, v, u )""",
'B,psi,u' :
"""dw_piezo_coupling.i1.Y2( inclusion.B, u, psi
)""",
'D,psi,phi' :
"""dw_diffusion.i1.Y2( inclusion.D, psi, phi )""",
'M,v,u' : """dw_mass_vector.i1.Y2( inclusion.density, v, u
)""",
}
...
mtx = assemble_by_blocks(conf.equations, problem,
ebcs = conf.ebcs,
epbcs = conf.epbcs)
Ah that's it! I did find out that the variables must appear in the
equation name, but I was setting my equations as actual equations
(i.e., """something = something""") and that wasn't
working properly
for some reason.
Is there a convention in SfePy as to how the equations are to be
interpreted when only their components are listed (as in your example
above), i.e., does SfePy assume that the first bit is a left-hand
side, the second bit is a right-hand side, and so on? Or is that left
to the programmer?
mtx is then a dict with 'K', ..., 'M'
keys. Of course, the corresponding
fields, materials, integrals, regions etc. have to be defined properly.
Thanks for all the help. Sorry for bombarding you
with questions.
You are welcome! I am glad this mailing list is alive (finally). Do not
hesitate to ask, the code is not by any means perfect, and I would like to know
about the various (often unexpected) ways people would like to use it.
Cheers,
Dominique
--
You received this message because you are subscribed to the Google Groups
"sfepy-devel" group.
To post to this group, send email to sfepy...(a)googlegroups.com.
To unsubscribe from this group, send email to sfepy-devel...(a)googlegroups.com.
For more options, visit this group at