Help expressing 1d eigenproblems in sfepy
Here's yet another plea for help on a "first problem".
I've been trying to learn FEM a little more thoroughly. Being a quantum chemist, I started thinking about looking at idealized 1d eigenproblems using FEA. I worked out some of the integrals, and the result is in the paper that I've added to the group space here:
http://groups.google.com/group/sfepy-devel/web/fem-1d-eigen.pdf
My quantum chemistry background shows through here, as I'm just using the finite elements like a normal basis set, a "pointy" Gaussian, if you will, and computing matrix elements. This isn't the way that real FEA people solve problems, and I'm trying to understand what that entails. I know it's something like:
- create a mesh
- work out a bunch of integrals over the mesh elements
- then a miracle occurs
- write down your final solution
Being a python person as well, I have high hopes that sfepy can help me understand how real FEA solutions work. But I'm still being a little bone-headed, and I don't understand the example files very well.
My hope is that the simple 1d problems that I've provided are simple enough that it wouldn't be too much work for someone (Ondrej?) to show me how to implement them in sfepy. Thanks in advance for any help anyone can offer.
Rick
I took a similar approach in working through a simple 1D problem and applying FEA to it. But it was fairly straight forward to implement mine in SfePy. I don't know what to say about yours, but once you get it figured out, I think it would make an excellent tutorial as well.
On Thu, Jul 31, 2008 at 3:33 PM, rpmu...@gmail.com <rpmu...@gmail.com> wrote:
Here's yet another plea for help on a "first problem".
I've been trying to learn FEM a little more thoroughly. Being a quantum chemist, I started thinking about looking at idealized 1d eigenproblems using FEA. I worked out some of the integrals, and the result is in the paper that I've added to the group space here:
http://groups.google.com/group/sfepy-devel/web/fem-1d-eigen.pdf
My quantum chemistry background shows through here, as I'm just using the finite elements like a normal basis set, a "pointy" Gaussian, if you will, and computing matrix elements. This isn't the way that real FEA people solve problems, and I'm trying to understand what that entails. I know it's something like:
- create a mesh
- work out a bunch of integrals over the mesh elements
- then a miracle occurs
- write down your final solution
Being a python person as well, I have high hopes that sfepy can help me understand how real FEA solutions work. But I'm still being a little bone-headed, and I don't understand the example files very well.
My hope is that the simple 1d problems that I've provided are simple enough that it wouldn't be too much work for someone (Ondrej?) to show me how to implement them in sfepy. Thanks in advance for any help anyone can offer.
Rick
On Thu, Jul 31, 2008 at 10:33 PM, rpmu...@gmail.com <rpmu...@gmail.com> wrote:
Here's yet another plea for help on a "first problem".
I've been trying to learn FEM a little more thoroughly. Being a quantum chemist, I started thinking about looking at idealized 1d eigenproblems using FEA. I worked out some of the integrals, and the result is in the paper that I've added to the group space here:
http://groups.google.com/group/sfepy-devel/web/fem-1d-eigen.pdf
My quantum chemistry background shows through here, as I'm just using the finite elements like a normal basis set, a "pointy" Gaussian, if you will, and computing matrix elements. This isn't the way that real FEA people solve problems, and I'm trying to understand what that entails. I know it's something like:
- create a mesh
- work out a bunch of integrals over the mesh elements
- then a miracle occurs
- write down your final solution
You do it exactly right in the paper. It's the same in 2D and 3D. I view FEM in the same way as you do, the mesh just defines a basis, you plug the basis in the equation, you get the matrix formulation, you solve it, you are done.
If you plugged the basis into the equation naively (e.g. directly), you'd get unsymmetric matrices, so you first integrate the equation per partes (or the green theorem, as you use in the paper).
Mathematicians also say a lot of other things about weak solutions and stuff, but I found it's not really useful and it doesn't provide any more useful information for me besides the above.
Being a python person as well, I have high hopes that sfepy can help me understand how real FEA solutions work. But I'm still being a little bone-headed, and I don't understand the example files very well.
My hope is that the simple 1d problems that I've provided are simple enough that it wouldn't be too much work for someone (Ondrej?) to show me how to implement them in sfepy. Thanks in advance for any help anyone can offer.
I am not sure if sfepy can solve 1D problems. 1D problems can be solved imho more efficiently by just integrating the ordinary differential equation. Ryan, do you know how to implement 1D stuff in sfepy?
For 2D and 3D, see the schroedinger.py script, that implements the particle in the box and oscillator problems, both in 2D and 3D. Maybe it could be adapted to 1D as well, I never thought about that.
Let me know if something is not clear in the script. If you are interested in the internals, I can also answer something that I understand and Robert can answer all the other details, as he programmed it in the first place.
Ondrej
Ondrej,
Thanks, I'll go through the schrodinger problems a little more carefully, and post my questions to this thread.
On Jul 31, 2:57 pm, "Ondrej Certik" <ond...@certik.cz> wrote:
On Thu, Jul 31, 2008 at 10:33 PM, rpmu...@gmail.com <rpmu...@gmail.com> wrote:
Here's yet another plea for help on a "first problem".
I've been trying to learn FEM a little more thoroughly. Being a quantum chemist, I started thinking about looking at idealized 1d eigenproblems using FEA. I worked out some of the integrals, and the result is in the paper that I've added to the group space here:
http://groups.google.com/group/sfepy-devel/web/fem-1d-eigen.pdf
My quantum chemistry background shows through here, as I'm just using the finite elements like a normal basis set, a "pointy" Gaussian, if you will, and computing matrix elements. This isn't the way that real FEA people solve problems, and I'm trying to understand what that entails. I know it's something like:
- create a mesh
- work out a bunch of integrals over the mesh elements
- then a miracle occurs
- write down your final solution
You do it exactly right in the paper. It's the same in 2D and 3D. I view FEM in the same way as you do, the mesh just defines a basis, you plug the basis in the equation, you get the matrix formulation, you solve it, you are done.
If you plugged the basis into the equation naively (e.g. directly), you'd get unsymmetric matrices, so you first integrate the equation per partes (or the green theorem, as you use in the paper).
Mathematicians also say a lot of other things about weak solutions and stuff, but I found it's not really useful and it doesn't provide any more useful information for me besides the above.
Being a python person as well, I have high hopes that sfepy can help me understand how real FEA solutions work. But I'm still being a little bone-headed, and I don't understand the example files very well.
My hope is that the simple 1d problems that I've provided are simple enough that it wouldn't be too much work for someone (Ondrej?) to show me how to implement them in sfepy. Thanks in advance for any help anyone can offer.
I am not sure if sfepy can solve 1D problems. 1D problems can be solved imho more efficiently by just integrating the ordinary differential equation. Ryan, do you know how to implement 1D stuff in sfepy?
For 2D and 3D, see the schroedinger.py script, that implements the particle in the box and oscillator problems, both in 2D and 3D. Maybe it could be adapted to 1D as well, I never thought about that.
Let me know if something is not clear in the script. If you are interested in the internals, I can also answer something that I understand and Robert can answer all the other details, as he programmed it in the first place.
Ondrej
I don't actually do 1D, but the other dimensions in my problem are doing Poisson-type displacement. So, only 1D is really active, but all are calculated.
On Thu, Jul 31, 2008 at 4:07 PM, rpmu...@gmail.com <rpmu...@gmail.com> wrote:
Ondrej,
Thanks, I'll go through the schrodinger problems a little more carefully, and post my questions to this thread.
On Jul 31, 2:57 pm, "Ondrej Certik" <ond...@certik.cz> wrote:
On Thu, Jul 31, 2008 at 10:33 PM, rpmu...@gmail.com <rpmu...@gmail.com> wrote:
Here's yet another plea for help on a "first problem".
I've been trying to learn FEM a little more thoroughly. Being a quantum chemist, I started thinking about looking at idealized 1d eigenproblems using FEA. I worked out some of the integrals, and the result is in the paper that I've added to the group space here:
http://groups.google.com/group/sfepy-devel/web/fem-1d-eigen.pdf
My quantum chemistry background shows through here, as I'm just using the finite elements like a normal basis set, a "pointy" Gaussian, if you will, and computing matrix elements. This isn't the way that real FEA people solve problems, and I'm trying to understand what that entails. I know it's something like:
- create a mesh
- work out a bunch of integrals over the mesh elements
- then a miracle occurs
- write down your final solution
You do it exactly right in the paper. It's the same in 2D and 3D. I view FEM in the same way as you do, the mesh just defines a basis, you plug the basis in the equation, you get the matrix formulation, you solve it, you are done.
If you plugged the basis into the equation naively (e.g. directly), you'd get unsymmetric matrices, so you first integrate the equation per partes (or the green theorem, as you use in the paper).
Mathematicians also say a lot of other things about weak solutions and stuff, but I found it's not really useful and it doesn't provide any more useful information for me besides the above.
Being a python person as well, I have high hopes that sfepy can help me understand how real FEA solutions work. But I'm still being a little bone-headed, and I don't understand the example files very well.
My hope is that the simple 1d problems that I've provided are simple enough that it wouldn't be too much work for someone (Ondrej?) to show me how to implement them in sfepy. Thanks in advance for any help anyone can offer.
I am not sure if sfepy can solve 1D problems. 1D problems can be solved imho more efficiently by just integrating the ordinary differential equation. Ryan, do you know how to implement 1D stuff in sfepy?
For 2D and 3D, see the schroedinger.py script, that implements the particle in the box and oscillator problems, both in 2D and 3D. Maybe it could be adapted to 1D as well, I never thought about that.
Let me know if something is not clear in the script. If you are interested in the internals, I can also answer something that I understand and Robert can answer all the other details, as he programmed it in the first place.
Ondrej
Hi Rick!
rpmu...@gmail.com wrote:
Here's yet another plea for help on a "first problem".
I've been trying to learn FEM a little more thoroughly. Being a quantum chemist, I started thinking about looking at idealized 1d eigenproblems using FEA. I worked out some of the integrals, and the result is in the paper that I've added to the group space here:
http://groups.google.com/group/sfepy-devel/web/fem-1d-eigen.pdf
My quantum chemistry background shows through here, as I'm just using the finite elements like a normal basis set, a "pointy" Gaussian, if you will, and computing matrix elements. This isn't the way that real FEA people solve problems, and I'm trying to understand what that entails. I know it's something like:
- create a mesh
- work out a bunch of integrals over the mesh elements
- then a miracle occurs
- write down your final solution
You got it completely right, it's just that. The miracle is in using those "pointy" Gaussians, or 'base functions with compact support' :) That way a global base function can be constructed as a sum of local, element-based, base functions - this is what the FE assembling is about.
Being a python person as well, I have high hopes that sfepy can help me understand how real FEA solutions work. But I'm still being a little bone-headed, and I don't understand the example files very well.
Look at input/poisson.py, that has some comments. Tell us, please, if they are not clear enough. We want them to be readable for FE newcomers like you.
My hope is that the simple 1d problems that I've provided are simple enough that it wouldn't be too much work for someone (Ondrej?) to show me how to implement them in sfepy. Thanks in advance for any help anyone can offer.
SfePy currently does not do 1D, just 2D or 3D (well, "just" is not the right word). It would not be difficult to extend it, though it is not probably worth it - standard ODE solvers are quite good.
Ondrej Certik wrote:
On Thu, Jul 31, 2008 at 10:33 PM, rpmu...@gmail.com <rpmu...@gmail.com> wrote:
http://groups.google.com/group/sfepy-devel/web/fem-1d-eigen.pdf
My quantum chemistry background shows through here, as I'm just using the finite elements like a normal basis set, a "pointy" Gaussian, if you will, and computing matrix elements. This isn't the way that real FEA people solve problems, and I'm trying to understand what that entails. I know it's something like:
- create a mesh
- work out a bunch of integrals over the mesh elements
- then a miracle occurs
- write down your final solution
You do it exactly right in the paper. It's the same in 2D and 3D. I view FEM in the same way as you do, the mesh just defines a basis, you plug the basis in the equation, you get the matrix formulation, you solve it, you are done.
If you plugged the basis into the equation naively (e.g. directly), you'd get unsymmetric matrices, so you first integrate the equation per partes (or the green theorem, as you use in the paper).
We do not care so much about symmetry (we use UMFPACK, right?), but about the necessary smoothness of the function spaces. Laplacian requires C^2 - now you cannot use P1 (linear) elements for that! This is why per-partes and the integral formulation is done.
Mathematicians also say a lot of other things about weak solutions and stuff, but I found it's not really useful and it doesn't provide any more useful information for me besides the above.
It does.
My hope is that the simple 1d problems that I've provided are simple enough that it wouldn't be too much work for someone (Ondrej?) to show me how to implement them in sfepy. Thanks in advance for any help anyone can offer.
I am not sure if sfepy can solve 1D problems. 1D problems can be solved imho more efficiently by just integrating the ordinary differential equation. Ryan, do you know how to implement 1D stuff in sfepy?
Right, it is not there, and is not probably a worthy addition.
For 2D and 3D, see the schroedinger.py script, that implements the particle in the box and oscillator problems, both in 2D and 3D. Maybe it could be adapted to 1D as well, I never thought about that.
It could, do you want to take a stab? :)
Let me know if something is not clear in the script. If you are interested in the internals, I can also answer something that I understand and Robert can answer all the other details, as he programmed it in the first place.
I will answer, but my answers tend to be somewhat time-lagged. I do not carry my laptop all the time like you do ;)
r.
If you plugged the basis into the equation naively (e.g. directly), you'd get unsymmetric matrices, so you first integrate the equation per partes (or the green theorem, as you use in the paper).
We do not care so much about symmetry (we use UMFPACK, right?), but about the necessary smoothness of the function spaces. Laplacian requires C^2 - now you cannot use P1 (linear) elements for that! This is why per-partes and the integral formulation is done.
Yes, that is the second reason for that.
Mathematicians usually say that the weak solution contains more solutions than the original equation, but obviously this is not true.
Ondrej
Ondrej Certik wrote:
If you plugged the basis into the equation naively (e.g. directly), you'd get unsymmetric matrices, so you first integrate the equation per partes (or the green theorem, as you use in the paper). We do not care so much about symmetry (we use UMFPACK, right?), but about the necessary smoothness of the function spaces. Laplacian requires C^2 - now you cannot use P1 (linear) elements for that! This is why per-partes and the integral formulation is done.
Yes, that is the second reason for that.
Mathematicians usually say that the weak solution contains more solutions than the original equation, but obviously this is not true.
H^1 is bigger than C^2, so in this sense, yes. They mean that the choice of the possible shapes is bigger. The actual equation, no matter the formulation should have just one solution to be well-posed, true.
r.
On Fri, Aug 1, 2008 at 12:23 PM, Ondrej Certik <ond...@certik.cz> wrote:
If you plugged the basis into the equation naively (e.g. directly), you'd get unsymmetric matrices, so you first integrate the equation per partes (or the green theorem, as you use in the paper).
We do not care so much about symmetry (we use UMFPACK, right?), but
BTW, we do actually care a lot, because in my case we are using an eigensolver and if the matrix is not symmetric, the eigenvalues are in general complex. Not mentioning that pysparse (imho) can only solve symmetric matrices. So using the per partes is actually crutial just to be able to solve it easily.
The argument about linear elements is imho not that important, you can always use higher order elements.
So from my point of view, you really want to rewrite the problem as symmetric just to move forward. Everything else is just a sideeffect, that can be fixed somehow.
Ondrej
Ondrej Certik wrote:
On Fri, Aug 1, 2008 at 12:23 PM, Ondrej Certik <ond...@certik.cz> wrote:
If you plugged the basis into the equation naively (e.g. directly), you'd get unsymmetric matrices, so you first integrate the equation per partes (or the green theorem, as you use in the paper). We do not care so much about symmetry (we use UMFPACK, right?), but
BTW, we do actually care a lot, because in my case we are using an eigensolver and if the matrix is not symmetric, the eigenvalues are in general complex. Not mentioning that pysparse (imho) can only solve symmetric matrices. So using the per partes is actually crutial just to be able to solve it easily.
Yes, in your case we do care.
The argument about linear elements is imho not that important, you can always use higher order elements.
But as you said, it would not be that easy.
So from my point of view, you really want to rewrite the problem as symmetric just to move forward. Everything else is just a sideeffect, that can be fixed somehow.
Symmetry is a win if one can achieve it, I agree with that. But believe me, weak formulations have other merits as well :)
r.
Robert Cimrman wrote:
Ondrej Certik wrote:
On Fri, Aug 1, 2008 at 12:23 PM, Ondrej Certik <ond...@certik.cz> wrote:
If you plugged the basis into the equation naively (e.g. directly), you'd get unsymmetric matrices, so you first integrate the equation per partes (or the green theorem, as you use in the paper). We do not care so much about symmetry (we use UMFPACK, right?), but BTW, we do actually care a lot, because in my case we are using an eigensolver and if the matrix is not symmetric, the eigenvalues are in general complex. Not mentioning that pysparse (imho) can only solve symmetric matrices. So using the per partes is actually crutial just to be able to solve it easily.
Yes, in your case we do care.
The argument about linear elements is imho not that important, you can always use higher order elements.
But as you said, it would not be that easy.
So from my point of view, you really want to rewrite the problem as symmetric just to move forward. Everything else is just a sideeffect, that can be fixed somehow.
Symmetry is a win if one can achieve it, I agree with that. But believe me, weak formulations have other merits as well :)
r.
Hi,
Let me chime in as a civil engineer. :-)
The 2nd order ODE for an axially loaded rod given by Ryan could also represent a tensioned string subject to lateral loads. For a point load, the deflected shape will be 2 straight line segments which has a discontinuity of slope at the location of the point load and hence not C^2. The weak formulation incorporates this very naturally. Further more, weak formulation is equivalent to the principle of virtual work in structural mechanics.
My 2 cents, ST
On Fri, Aug 1, 2008 at 4:54 PM, LUK ShunTim <shunt...@polyu.edu.hk> wrote:
Robert Cimrman wrote:
Ondrej Certik wrote:
On Fri, Aug 1, 2008 at 12:23 PM, Ondrej Certik <ond...@certik.cz> wrote:
If you plugged the basis into the equation naively (e.g. directly), you'd get unsymmetric matrices, so you first integrate the equation per partes (or the green theorem, as you use in the paper). We do not care so much about symmetry (we use UMFPACK, right?), but BTW, we do actually care a lot, because in my case we are using an eigensolver and if the matrix is not symmetric, the eigenvalues are in general complex. Not mentioning that pysparse (imho) can only solve symmetric matrices. So using the per partes is actually crutial just to be able to solve it easily.
Yes, in your case we do care.
The argument about linear elements is imho not that important, you can always use higher order elements.
But as you said, it would not be that easy.
So from my point of view, you really want to rewrite the problem as symmetric just to move forward. Everything else is just a sideeffect, that can be fixed somehow.
Symmetry is a win if one can achieve it, I agree with that. But believe me, weak formulations have other merits as well :)
r.
Hi,
Let me chime in as a civil engineer. :-)
The 2nd order ODE for an axially loaded rod given by Ryan could also represent a tensioned string subject to lateral loads. For a point load, the deflected shape will be 2 straight line segments which has a discontinuity of slope at the location of the point load and hence not C^2. The weak formulation incorporates this very naturally. Further more, weak formulation is equivalent to the principle of virtual work in structural mechanics.
I don't know the equations in this case, but I guess it's the same as in electromagnetism, where you also get solutions that are discontinuous at the boundary (either tangent or perpedicular projection) and no one doubts that maxwells equations govern it. I mean I don't understand why you require C^2? The discontinuity in slope means that the derivative is not continuous, right? E.g. the laplace operator generates couple delta functions, correct? I have no problems with that and I don't need any weak formulation for that. C^2 is just something that mathematicians invented, so that they have some justifications to think that weak formulation is something more, but in fact everybody knows that the solution that you are looking for is not C^2 and you can get it from the original differential equation as easily, or is there something more that I don't see?
Ondrej
Thanks for the reply, Robert.
The reason I think that I'm doing something wrong is that I worked out all of my integrals in terms of the element functions themselves, whereas most of the papers I've read (for example, ) work out the integrals in terms of shape functions on [0,1], typically using Gaussian quadrature, and then assemble the matrix elements over the finite elements.
On Aug 1, 4:11 am, Robert Cimrman <cimr...@ntc.zcu.cz> wrote:
Hi Rick!
rpmu...@gmail.com wrote:
Here's yet another plea for help on a "first problem".
I've been trying to learn FEM a little more thoroughly. Being a quantum chemist, I started thinking about looking at idealized 1d eigenproblems using FEA. I worked out some of the integrals, and the result is in the paper that I've added to the group space here:
http://groups.google.com/group/sfepy-devel/web/fem-1d-eigen.pdf
My quantum chemistry background shows through here, as I'm just using the finite elements like a normal basis set, a "pointy" Gaussian, if you will, and computing matrix elements. This isn't the way that real FEA people solve problems, and I'm trying to understand what that entails. I know it's something like:
- create a mesh
- work out a bunch of integrals over the mesh elements
- then a miracle occurs
- write down your final solution
You got it completely right, it's just that. The miracle is in using those "pointy" Gaussians, or 'base functions with compact support' :) That way a global base function can be constructed as a sum of local, element-based, base functions - this is what the FE assembling is about.
Being a python person as well, I have high hopes that sfepy can help me understand how real FEA solutions work. But I'm still being a little bone-headed, and I don't understand the example files very well.
Look at input/poisson.py, that has some comments. Tell us, please, if they are not clear enough. We want them to be readable for FE newcomers like you.
My hope is that the simple 1d problems that I've provided are simple enough that it wouldn't be too much work for someone (Ondrej?) to show me how to implement them in sfepy. Thanks in advance for any help anyone can offer.
SfePy currently does not do 1D, just 2D or 3D (well, "just" is not the right word). It would not be difficult to extend it, though it is not probably worth it - standard ODE solvers are quite good.
I had meant to cite a paper below in the parentheses "(for example, )", but I couldn't find a version of it online. Hopefully people know what I'm talking about, though.
On Aug 1, 1:26 pm, "rpmu...@gmail.com" <rpmu...@gmail.com> wrote:
Thanks for the reply, Robert.
The reason I think that I'm doing something wrong is that I worked out all of my integrals in terms of the element functions themselves, whereas most of the papers I've read (for example, ) work out the integrals in terms of shape functions on [0,1], typically using Gaussian quadrature, and then assemble the matrix elements over the finite elements.
On Aug 1, 4:11 am, Robert Cimrman <cimr...@ntc.zcu.cz> wrote:
Hi Rick!
rpmu...@gmail.com wrote:
Here's yet another plea for help on a "first problem".
I've been trying to learn FEM a little more thoroughly. Being a quantum chemist, I started thinking about looking at idealized 1d eigenproblems using FEA. I worked out some of the integrals, and the result is in the paper that I've added to the group space here:
http://groups.google.com/group/sfepy-devel/web/fem-1d-eigen.pdf
My quantum chemistry background shows through here, as I'm just using the finite elements like a normal basis set, a "pointy" Gaussian, if you will, and computing matrix elements. This isn't the way that real FEA people solve problems, and I'm trying to understand what that entails. I know it's something like:
- create a mesh
- work out a bunch of integrals over the mesh elements
- then a miracle occurs
- write down your final solution
You got it completely right, it's just that. The miracle is in using those "pointy" Gaussians, or 'base functions with compact support' :) That way a global base function can be constructed as a sum of local, element-based, base functions - this is what the FE assembling is about.
Being a python person as well, I have high hopes that sfepy can help me understand how real FEA solutions work. But I'm still being a little bone-headed, and I don't understand the example files very well.
Look at input/poisson.py, that has some comments. Tell us, please, if they are not clear enough. We want them to be readable for FE newcomers like you.
My hope is that the simple 1d problems that I've provided are simple enough that it wouldn't be too much work for someone (Ondrej?) to show me how to implement them in sfepy. Thanks in advance for any help anyone can offer.
SfePy currently does not do 1D, just 2D or 3D (well, "just" is not the right word). It would not be difficult to extend it, though it is not probably worth it - standard ODE solvers are quite good.
On Fri, Aug 1, 2008 at 9:26 PM, rpmu...@gmail.com <rpmu...@gmail.com> wrote:
Thanks for the reply, Robert.
The reason I think that I'm doing something wrong is that I worked out all of my integrals in terms of the element functions themselves, whereas most of the papers I've read (for example, ) work out the integrals in terms of shape functions on [0,1], typically using Gaussian quadrature, and then assemble the matrix elements over the finite elements.
Yes. It's the same thing though. Very brief explanation is in the draft of my master thesis I sent you. Instead of calculating the integral of basis functions directly, you first convert it to an integral over a reference element (introducing the Jacobian in there) and then to a sum over gauss points. This can be then easily automated.
Ondrej
Ondrej Certik wrote:
On Fri, Aug 1, 2008 at 9:26 PM, rpmu...@gmail.com <rpmu...@gmail.com> wrote:
Thanks for the reply, Robert.
The reason I think that I'm doing something wrong is that I worked out all of my integrals in terms of the element functions themselves, whereas most of the papers I've read (for example, ) work out the integrals in terms of shape functions on [0,1], typically using Gaussian quadrature, and then assemble the matrix elements over the finite elements.
Yes. It's the same thing though. Very brief explanation is in the draft of my master thesis I sent you. Instead of calculating the integral of basis functions directly, you first convert it to an integral over a reference element (introducing the Jacobian in there) and then to a sum over gauss points. This can be then easily automated.
Ondrej
Correct. Look also for 'isoparametric elements' on google: -> http://www.mathsoft.cse.clrc.ac.uk/felib/Docs/html/Intro/intro-node31.html ...
r.
Ondrej Certik wrote:
On Fri, Aug 1, 2008 at 4:54 PM, LUK ShunTim <shunt...@polyu.edu.hk> wrote:
Robert Cimrman wrote:
Ondrej Certik wrote:
On Fri, Aug 1, 2008 at 12:23 PM, Ondrej Certik <ond...@certik.cz> wrote:
> If you plugged the basis into the equation naively (e.g. directly), > you'd get unsymmetric matrices, so you first integrate the equation > per partes (or the green theorem, as you use in the paper). We do not care so much about symmetry (we use UMFPACK, right?), but BTW, we do actually care a lot, because in my case we are using an eigensolver and if the matrix is not symmetric, the eigenvalues are in general complex. Not mentioning that pysparse (imho) can only solve symmetric matrices. So using the per partes is actually crutial just to be able to solve it easily. Yes, in your case we do care.
The argument about linear elements is imho not that important, you can always use higher order elements. But as you said, it would not be that easy.
So from my point of view, you really want to rewrite the problem as symmetric just to move forward. Everything else is just a sideeffect, that can be fixed somehow. Symmetry is a win if one can achieve it, I agree with that. But believe me, weak formulations have other merits as well :)
r. Hi,
Let me chime in as a civil engineer. :-)
The 2nd order ODE for an axially loaded rod given by Ryan could also represent a tensioned string subject to lateral loads. For a point load, the deflected shape will be 2 straight line segments which has a discontinuity of slope at the location of the point load and hence not C^2. The weak formulation incorporates this very naturally. Further more, weak formulation is equivalent to the principle of virtual work in structural mechanics.
I don't know the equations in this case, but I guess it's the same as in electromagnetism, where you also get solutions that are discontinuous at the boundary (either tangent or perpedicular projection) and no one doubts that maxwells equations govern it. I mean I don't understand why you require C^2? The discontinuity in slope means that the derivative is not continuous, right? E.g. the laplace operator generates couple delta functions, correct? I have no problems with that and I don't need any weak formulation for that. C^2 is just something that mathematicians invented, so that they have some justifications to think that weak formulation is something more, but in fact everybody knows that the solution that you are looking for is not C^2 and you can get it from the original differential equation as easily, or is there something more that I don't see?
What is the definition of the (second) derivative you use? Once you use something with integrals and test functions (e.g. distributions), you have the weak formulation. You just do not like the term "weak", right? Well, it's just a definition of the exactly same thing you think about the solution.
Classical solution = C^2 Weak solution = H^1
It's a naming convention, and everybody uses it ;)
r.
participants (5)
-
LUK ShunTim -
Ondrej Certik -
Robert Cimrman -
rpmu...＠gmail.com -
Ryan Krauss