femlisp-user
[Top][All Lists]

## Re: [femlisp-user] Constructing 2D domain

 From: Nicolas Neuss Subject: Re: [femlisp-user] Constructing 2D domain Date: Sun, 07 Nov 2010 17:59:24 +0100 User-agent: Gnus/5.13 (Gnus v5.13) Emacs/23.2 (gnu/linux)

```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)

```