Re: [femlisp-user] Constructing 2D domain

[Sebastian Sturm]

First of all, thank you very much for your detailed solution! Looking at the 
source code, I would definitely not have been able to achieve this by guessing. 
My apologies for replying so late; the reason is that I'm encountering a 
probably embarrassingly trivial problem that, up to now, prevents me from 
looking at the solution. When using (plot :program :vtk) to export the 
solution, SBCL crashes inside fl.graphic::graphic-write-data. For the 
ellsys-model-problem, (plot :program :vtk) works as expected, though. I thought 
I would be able to access the solution manually by calling something like 
(fe-value (getbb *result* :solution) (vector 0.5 0.5)), but this produces the 
error message
'There is no applicable method for the generic function 
#<standard-generic-function fl.mesh:global->local (3)> when called with 
arguments (nil #(0.5 0.5))'. Since the test has passed on your machine, 
probably something is wrong with my femlisp installation. I'll try to figure 
out what's broken, just wanted to let you know why I'm taking so long to answer.

thanks again,
Sebastian
On 07.11.2010, at 17:59, Nicolas Neuss wrote:

> Sebastian Sturm <address@hidden> writes:
> 
>>> For this application, I used something like the following:
>>> 
>>> <sebastian-sturm.lisp>
>> 
>> Great, thank you very much! That is exactly what I have been searching
>> for. I now have two separate domains (outer and inner part) that
>> should have different material parameters sigma_1 and sigma_2,
>> respectively. I would like to solve div(sigma(x) * grad(phi)) = 0 on
>> the sum of both domains, with phi(x) at the inner (circular hole)
>> boundaries set to +1 and -1, respectively, phi(x) = 0 or grad(phi(x))
>> = 0 at the outermost boundary and sigma_1 (d phi/dn)_{boundary of
>> domain 1} = sigma_2 (d phi/dn)_{boundary of domain 2} at the interface
>> between the inner domain and the outer domain. Is that possible?
>> 
>> I have attached a pdf file showing the domain I'm referring to. From
>> browsing the source code, I guess I should set up an elliptic system,
>> define ellsys::a as a simple identity matrix with diagonal elements
>> sigma(x) and (somehow) set up Dirichlet boundary conditions at the
>> holes. I assume I can figure that out on my own, but identifying the
>> inner boundary of the outer domain with the outer boundary of the
>> inner domain and setting up the right BC seems more daunting to me.
> 
> Probably you had the right feeling.  It is not too difficult, but the
> setup is also not completely trivial.  At the moment, Femlisp does not
> have a GUI like COMSOL, for example, where you can graphically draw
> regions and specify boundary conditions.  I could surely add a special
> DSL for such stuff.  But it gets difficult for 3D domains, systems of
> PDEs, and I have the feeling that I still need some experience before
> getting this completely right.
> 
>> I also guess I could do without the BC at the interface region by
>> constructing a finite-width boundary region using telescope and
>> explicitly interpolating sigma(x) between its two extreme values
>> sigma_1 and sigma_2, in case the current ellsys model is not well
>> suited to the boundary conditions described above. In principle,
>> however, I would very much like to know how to set up the most general
>> kind of interface/boundary conditions.
> 
> For the boundary conditions you described above, you don't have to do
> anything for a finite element discretizations, because they are
> "natural", i.e. the FE solution satisfies them automatically.
> 
>> My apologies for stealing your time, but my weak Lisp skills do not
>> cater very well to the trial-and-error approach. In compensation, I
>> could try to contribute a tutorial to the femlisp manual explaining
>> these basic techniques (once I have actually understood them) to
>> newcomers like me.
> 
> If you feel like it, this would probably be a good addition.  I have
> programmed a Femlisp demo which does mostly what you want.  As much as I
> see, the only difference is that I specified the outer and inner
> boundary as squares for simplicity (the electrodes are circles in the
> innermost square however).  It should be completely straightforward to
> change this to your H-shape (and I would be interested in the precise
> data for making the application more realistic).  Maybe you'll also have
> to change some parameters for the mesh generation.
> 
> You can get that code by doing a "cvs update" in the Femlisp directory
> again.  After that you can run the demo by doing (femlisp-demo) and then
> choosing "Equations->Laplace->Electromagnetic Potential".
> 
> You'll find the source code in
> #p"femlisp:src;applications;cdr;sturm.lisp".
> 
> Nicolas
> 
> P.S.: Note that the code is not highly performant, however.  To achieve
> that would be a more elaborate task consisting of
> 
> (a) generating a better mesh (probably possible only by avoiding
>   Triangle)
> 
> (b) choosing an adaptive algorithm with a suitable error estimator (for
>   this one needs to know especially the quantity one is interested in)
> 
> _______________________________________________
> femlisp-user mailing list
> address@hidden
> http://lists.nongnu.org/mailman/listinfo/femlisp-user