PDE toolbox
Ondrej Certik
On Dec 7, 2007 6:50 AM, Joshua Kantor <kant...@gmail.com> wrote:
>
> I would be interested in helping with a PDE toolbox. I didn't want to
> work on it alone as I'm pretty sure I'd make some stupid design
> choices. It would be nice to start some work on PDE functionality in
> SAGE.
I would be very interested in this. Especially with finite elements
solver. SciPy has some solvers, but even producing the element
matrices would
be enough.
See also:
http://groups.google.com/group/sympy/browse_thread/thread/f73b00d1bdf0255a
http://groups.google.com/group/sympy/browse_thread/thread/6b8b4cbaa1914c45
and the project SymFe:
http://code.google.com/p/symfe/
Basically, the way the equation is defined just need SAGE.calculus.
But I would very much like this to be modular, so that we can use both
maxima and sympy as a backend, so that I can use this in my own
programs, because depending on the whole SAGE would be too much.
However, it's a huge task. There are many related projects, mainly:
http://www.fenics.org/wiki/FEniCS_Project
http://www.fenics.org/wiki/SyFi
But let's discuss at least the design. Let's use the same design as in:
http://www.mathworks.com/products/pde/
or is there a better product out there?
Ondrej
Tim Lahey
On Dec 7, 2007, at 4:26 AM, Ondrej Certik wrote:
>
> Hi,
>
> On Dec 7, 2007 6:50 AM, Joshua Kantor <kant...@gmail.com> wrote:
>>
>> I would be interested in helping with a PDE toolbox. I didn't want to
>> work on it alone as I'm pretty sure I'd make some stupid design
>> choices. It would be nice to start some work on PDE functionality in
>> SAGE.
>
> I would be very interested in this. Especially with finite elements
> solver. SciPy has some solvers, but even producing the element
> matrices would
> be enough.
>
For variational problems, I've already written code in Maple to derive
the element matrices. I presented a paper on in at the Maple 2005
Conference. Unfortunately, I've been having a difficult time translating
some of the things I did from Maple to Sage. It was only for 1D problems
(since that was all I needed), but it is modular enough to extend
to others. The code derives the mass and stiffness matrices and then
outputs them as MATLAB functions. I wrote the code because I needed
something that can handle rotating flexible beams.
>
>
> Basically, the way the equation is defined just need SAGE.calculus.
> But I would very much like this to be modular, so that we can use both
> maxima and sympy as a backend, so that I can use this in my own
> programs, because depending on the whole SAGE would be too much.
>
> However, it's a huge task. There are many related projects, mainly:
>
> http://www.fenics.org/wiki/FEniCS_Project
> http://www.fenics.org/wiki/SyFi
I know about both of these, but I had never heard of the symfe work.
Interesting, since my Maple toolbox is called SFEM.
There is also a package for Maple called femLego:
http://www.mech.kth.se/~gustava/femLego/
It's design is about solving the PDE once it has been derived. It
uses Maple to write Fortran code, and compiles it with its own code.
I nearly ended up using this, but see below.
> But let's discuss at least the design. Let's use the same design as
> in:
>
> http://www.mathworks.com/products/pde/
There are some problems with its design. It is basically an early
version
of FEMLAB, and it puts its focus strictly on solving PDEs once the
formulation of the PDE has already been done. If a Ritz approach is
used,
the formulation of the PDE is never explicitly done and this design
can't
handle it. Plus, for systems with multiple coupled variables, the
Galerkin
approach becomes very complicated. I know because my SFEM package was
originally an attempt to derive the weak form for use in FEMLAB.
>
> or is there a better product out there?
>
> Ondrej
The difficulty is that most products seem to be trying to replace the
major
FE packages, but even those who do symbolic work, don't put any focus on
the derivation. If one is going to use the weak formulation, this can be
difficult.
Tim.
Ondrej Certik
> the element matrices. I presented a paper on in at the Maple 2005
> Conference. Unfortunately, I've been having a difficult time translating
> some of the things I did from Maple to Sage. It was only for 1D problems
> (since that was all I needed), but it is modular enough to extend
> to others. The code derives the mass and stiffness matrices and then
> outputs them as MATLAB functions. I wrote the code because I needed
> something that can handle rotating flexible beams.
Is your matlab code available? Where did you have problems translating
it to Python?
> I know about both of these, but I had never heard of the symfe work.
> Interesting, since my Maple toolbox is called SFEM.
That's because we started symfe a week ago. :) Just experimenting.
> There is also a package for Maple called femLego:
>
> http://www.mech.kth.se/~gustava/femLego/
>
> It's design is about solving the PDE once it has been derived. It
> uses Maple to write Fortran code, and compiles it with its own code.
> I nearly ended up using this, but see below.
>
> > But let's discuss at least the design. Let's use the same design as
> > in:
> >
> > http://www.mathworks.com/products/pde/
>
> There are some problems with its design. It is basically an early
> version
Are you talking about the femLego, or pde from mathworks?
> of FEMLAB, and it puts its focus strictly on solving PDEs once the
> formulation of the PDE has already been done. If a Ritz approach is
> used,
> the formulation of the PDE is never explicitly done and this design
> can't
> handle it. Plus, for systems with multiple coupled variables, the
> Galerkin
> approach becomes very complicated. I know because my SFEM package was
> originally an attempt to derive the weak form for use in FEMLAB.
Where is the SFEM package? Is it this:
http://people.civil.gla.ac.uk/~bordas/sfem.html
?
> > or is there a better product out there?
> >
> > Ondrej
>
> The difficulty is that most products seem to be trying to replace the
> major
> FE packages, but even those who do symbolic work, don't put any focus on
> the derivation. If one is going to use the weak formulation, this can be
> difficult.
I don't understand what you mean by "the derivation" - like taking the
differential equation and constructing the corresponding integral
(variational) formulation?
Ondrej
Tim Lahey
On Dec 8, 2007, at 8:32 AM, Ondrej Certik wrote:
>
> Is your matlab code available? Where did you have problems translating
> it to Python?
>
My MATLAB code isn't available at the moment. It shouldn't be a
problem translating it to Python. The problem is translating the Maple
code that derives the element matrices to Sage. I'm not exactly thrilled
with the design of the MATLAB code, but it works for my problems. I'd
probably make it object-oriented when converting it to Python.
>> I know about both of these, but I had never heard of the symfe work.
>> Interesting, since my Maple toolbox is called SFEM.
>
> That's because we started symfe a week ago. :) Just experimenting.
Ah, that's interesting. I thought it was just because I hadn't looked in
a while.
>
>> There is also a package for Maple called femLego:
>>
>> http://www.mech.kth.se/~gustava/femLego/
>>
>> It's design is about solving the PDE once it has been derived. It
>> uses Maple to write Fortran code, and compiles it with its own code.
>> I nearly ended up using this, but see below.
>>
>>> But let's discuss at least the design. Let's use the same design as
>>> in:
>>>
>>> http://www.mathworks.com/products/pde/
>>
>> There are some problems with its design. It is basically an early
>> version
>
> Are you talking about the femLego, or pde from mathworks?
>
I'm talking about the PDE toolbox. I don't have too many problems with
femLego other than it assumes you already have the weak form.
>> of FEMLAB, and it puts its focus strictly on solving PDEs once the
>> formulation of the PDE has already been done. If a Ritz approach is
>> used,
>> the formulation of the PDE is never explicitly done and this design
>> can't
>> handle it. Plus, for systems with multiple coupled variables, the
>> Galerkin
>> approach becomes very complicated. I know because my SFEM package was
>> originally an attempt to derive the weak form for use in FEMLAB.
>
> Where is the SFEM package? Is it this:
>
> http://people.civil.gla.ac.uk/~bordas/sfem.html
>
No, the SFEM package I wrote isn't currently available. I was planning
on waiting on releasing the code (both the Maple and MATLAB portions)
until after I finished a paper on the parts for the Journal of Symbolic
Computation. I didn't realize that somebody else had a package out there
with that name.
>
>>> or is there a better product out there?
>>>
>>> Ondrej
>>
>> The difficulty is that most products seem to be trying to replace the
>> major
>> FE packages, but even those who do symbolic work, don't put any
>> focus on
>> the derivation. If one is going to use the weak formulation, this
>> can be
>> difficult.
>
> I don't understand what you mean by "the derivation" - like taking the
> differential equation and constructing the corresponding integral
> (variational) formulation?
>
What I mean by the derivation is the steps to get the finite element
approximation. There are two basic approaches to this.
1. Galerkin method - Convert the PDE to the weak form.
2. Ritz method - Requires a variational statement for the underlying
physical problem (e.g., Hamilton's Principle). The PDE is never
explicitly
formed since the FE approximation is used directly in the variational
statement.
The Ritz method is less general, since you need an underlying
variational
principle, but it is easier to handle coupled equations. Getting the
weak form of the equations is difficult if the equations are coupled.
None of the introductory FE textbooks that I've seen covers it. Examples
always seem to use only one variable or multiple decoupled variables.
For my problems, the Ritz method is applicable so I use it.
Unfortunately,
most packages assume you use Galerkin. In addition, they assume you
already
have the weak equations and like I've stated, getting them can be
challenging.
So, given that this is part of a symbolic package, there should be some
support for that.
Tim.
Ondrej Certik
> problem translating it to Python. The problem is translating the Maple
> code that derives the element matrices to Sage. I'm not exactly thrilled
> with the design of the MATLAB code, but it works for my problems. I'd
> probably make it object-oriented when converting it to Python.
And which exact symbolic features are still missing in Sage, if any? I
mean, was the problem just that python is
different than Matlab, or just because something in Sage should be improved?
> What I mean by the derivation is the steps to get the finite element
> approximation. There are two basic approaches to this.
>
> 1. Galerkin method - Convert the PDE to the weak form.
>
> 2. Ritz method - Requires a variational statement for the underlying
> physical problem (e.g., Hamilton's Principle). The PDE is never
> explicitly
> formed since the FE approximation is used directly in the variational
> statement.
>
> The Ritz method is less general, since you need an underlying
> variational
> principle, but it is easier to handle coupled equations. Getting the
> weak form of the equations is difficult if the equations are coupled.
> None of the introductory FE textbooks that I've seen covers it. Examples
> always seem to use only one variable or multiple decoupled variables.
>
> For my problems, the Ritz method is applicable so I use it.
> Unfortunately,
> most packages assume you use Galerkin. In addition, they assume you
> already
> have the weak equations and like I've stated, getting them can be
> challenging.
> So, given that this is part of a symbolic package, there should be some
> support for that.
I got it now. Yes, I only need the Galerkin method. I have the
equation, I have the weak form. Now I want to solve it, on a general
geometry with general boundary conditions. I use libmesh for that:
http://libmesh.sourceforge.net/
it's a very good FEM library. But still maybe it could be done
simpler, using the SyFi approach or similar.
Ondrej
Tim Lahey
On Dec 8, 2007, at 9:13 AM, Ondrej Certik wrote:
>
>> My MATLAB code isn't available at the moment. It shouldn't be a
>> problem translating it to Python. The problem is translating the
>> Maple
>> code that derives the element matrices to Sage. I'm not exactly
>> thrilled
>> with the design of the MATLAB code, but it works for my problems. I'd
>> probably make it object-oriented when converting it to Python.
>
> And which exact symbolic features are still missing in Sage, if any? I
> mean, was the problem just that python is
> different than Matlab, or just because something in Sage should be
> improved?
I make heavy use of Maple's op command and hastype (I find the op
sequences
of the integrals and integrands). The code is pretty messy so I'll
refrain
from posting it here at the moment. I also make use of the fact you can
create new variable names by appending to old ones in my Euler-
Lagrange code.
This code is much faster than Maple's code. My Euler-Lagrange code is
as follows. In Maple this executes about 2 orders of magnitude faster
than
the built-in one (0.1 seconds vs. 10 seconds) and it preserves the
order.
I assume that time is the independent variable, but modifying that is
easy.
The conversion to new variables is because Maple can't take
derivatives with
respect to variables, only symbols.
EulerLagrange := proc(Lagrangian::anything, variables::list)
local num_list, qv_name, vel_var, qv_subs, qv_unsubs, Lagrange_subs1,
Lagrange_subs2, dL_dqv1, dL_dqv2, dL_dqv, dL_dqvt, dL_dq, dL_dq1,
dL_dq2, dL_dq3, q_name, q_subs, q_unsubs:
# create a list of indices from 1 to the number of variables
# used in the formulation
num_list := [seq(i,i=1..nops(variables))]:
# Define a list of generalized velocity and position variables
qv_name := map2(cat,qv,num_list):
q_name := map2(cat,q,num_list):
# Equate the time derivatives of the system variable to the
# generalized velocities and also define the reverse mapping
vel_var := map(diff,variables,t):
qv_subs := zip(equate,vel_var,qv_name):
qv_unsubs := zip(equate,qv_name,vel_var):
# Equate the generalized positions to the system variables
# and define the reverse mapping
q_subs := zip(equate,variables,q_name):
q_unsubs := zip(equate,q_name,variables):
# Convert the Lagrangian to the generalized position and velocity
variables
Lagrange_subs1 := subs(qv_subs,Lagrangian):
Lagrange_subs2 := subs(q_subs,Lagrange_subs1):
# Differentiate the Lagrangian with respect to the
# generalized velocities and positions
dL_dqv1 := map2(diff,Lagrange_subs2,qv_name):
dL_dq1 := map2(diff,Lagrange_subs2,q_name):
# Revert back to the system variables
dL_dq2 := map2(subs,qv_unsubs,dL_dq1):
dL_dqv2 := map2(subs,qv_unsubs,dL_dqv1):
dL_dqv := map2(subs,q_unsubs,dL_dqv2):
dL_dq := map2(subs,q_unsubs,dL_dq2):
dL_dqvt := map(diff,dL_dqv,t):
# Return the two components of the Euler-Lagrange Equation
return (dL_dqvt, dL_dq):
end proc:
The equate function used in zip is defined as:
equate := (x,y)->x=y:
I used it over `=` because it is slightly faster since Maple doesn't
have
to convert equate to a function.
I'm happy to work with someone to convert things over to Sage, but I
can't
seem to find a number of things I'm used to in Maple.
Tim.
Ondrej Certik
Awesome, nice application. Could you please post here some examples of
usage? Correct input and output, for three or four cases, so that I
can easily see what exactly the input and output should look like.
I'll then think how to port this to Sage. Note: I'll be on mountains
the next week, so I might do it later.
Ondrej