Subject: [FEniCS-users] Defining a Robin type boundary condition with Laplace equation

From: Thomas Ward <[email-masked]>
Date: Tue, 22 Sep 2009 16:35:17 +0100

Hello Fenics people,

this is my first post to this list , and I have just started to try to use
ffc and dolfin.
I am impressed with the concept and implementation of your project, thanks
you for making it available.

1.
I am attempting to solve a simple incompressible potential flow problem
(laplace equation) on a regular cube
with
 dphi/dx=omega at xmin
 dphi/dx=-omega at xmax
 dphi/dy=0 at ymin and ymax and zmin
dphi/dz=-omega*omega*phi/g at zmax

2.
I think I have grasped that to define xmin and xmax neuman conditions I need
something like

element = FiniteElement("Lagrange", "tetrahedron", 2)
v = TestFunction(element)
u = TrialFunction(element)
f = Function(element)
g = Function(element)

a = dot(grad(v), grad(u))*dx
L = v*f*dx+v*g*ds

and then enforce the neuman condtions in g.eval so that
int_over_xmin(v*omega)
int_over_xmax(-v*omega)

this is the code I am suggesting

        real eval ( const real* x ) const {
            if ( x[0] >1- DOLFIN_EPS )  return omega;
            else
              if ( x[0] < DOLFIN_EPS )  return -omega;
              else return 0;
        }

I am slightly worried by this however because I thought that the dphi/dn
boundary conditions will be enforced naturally by the variational form and
it seems like I am enforcing them again with this g.eval by the
       else return 0;
statement ??

3.
I can't see how to imlpenent the zmax boundary condition in dolfin, I think
the  weak form should be
    int_over_zmax(alpha * u *v )
       where alpha is a constant = omega * omega /g
but I can't work out how to fit this into ffc or dolfin.

4.
On a separate issue are there any infinite elements implemented in fenics,
if no then are there any plans to do so and how difficult would it be to do
so?

thanks for any help

Tom Ward.

---

From: Anders Logg <[email-masked]>
Date: Thu, 24 Sep 2009 21:48:19 +0200

For Robin type boundary conditions in Laplace, just integrate by parts
and replace -du/dn with whatever it should be:

 (v, -du/dn)_\partial\Omega = (v, g + gamma (u - u0))_\partial\Omega

if the Robin condition is

  -du/dn = g + gamma(u - u0)

For a small value of gamma, this reduces to a Neumann condition
-du/dn = g and for a large value, it reduces to a weakly enforced
Dirichlet condition u = u0.

You will get an integral involving the unknown u on the boundary. Just
keep that on the left-hand side and throw everything else on the
right-hand side.

Regarding infinite elements, I don't think we have any such plans.

--
Anders

On Tue, Sep 22, 2009 at 04:35:17PM +0100, Thomas Ward wrote:
> Hello Fenics people,
>
> this is my first post to this list , and I have just started to try to use ffc
> and dolfin.
> I am impressed with the concept and implementation of your project, thanks you
> for making it available.
>
>
> 1.
> boundary conditions will be enforced naturally by the variational form and it
> seems like I am enforcing them again with this g.eval by the
>        else return 0;
> statement ??
>
>
> 3.
> I can't see how to imlpenent the zmax boundary condition in dolfin, I think the
>   weak form should be
> On a separate issue are there any infinite elements implemented in fenics, if
> no then are there any plans to do so and how difficult would it be to do so?
>
>
> thanks for any help
>
> Tom Ward.

> _______________________________________________
> FEniCS-users mailing list
> [email-masked]
> http://fenics.org/mailman/listinfo/fenics-users

---

From: Thomas Ward <[email-masked]>
Date: Fri, 25 Sep 2009 06:58:15 +0100

Many Thanks Anders,

I am still slightly unclear about the free surface BC
  dphi/dz = -constant *phi at zmax
because in

    void Function::eval
<http://www.fenics.org/pub/documents/dolfin/dolfin-progr-reference/d2/d42/classdolfin_1_1Function.html#44fbf0b7b6bf28abb54e6d02f7bf54e3>(double*
values, const double* x, const ufc::cell& ufc_cell,uint
<http://www.fenics.org/pub/documents/dolfin/dolfin-progr-reference/df/de3/namespacedolfin.html#38ba062a16feed21dca09d99a2eed784>
cell_index) const;

I believe that "values" argument is the test function?
I think that after integration by parts the free surface BC becomes
  integral_over_surface_of_zmax(constant * trial_function *test_function )
?

i.e. I need to integrate with the trial function as well as the test
function for this free surface BC?

Anyway I have just found that the users used to post to the dofin-dev list
and there are some similar question to mine there which may be helpful.
I'll have a dig around in the archives of the dolfin-dev and see if I can
work it out,

regards

Tom

On Thu, Sep 24, 2009 at 8:48 PM, Anders Logg <[email-masked]> wrote:

> -----BEGIN PGP SIGNATURE-----
> Version: GnuPG v1.4.9 (GNU/Linux)
>
> iEYEARECAAYFAkq7zQMACgkQTuwUCDsYZdE8BACeJY1ltvoeWRDh26GmH4JOjiHa
> qLsAn0nLl2hnUxgWHqQ0Mk/MWAazTJeS
> =U4PD
> -----END PGP SIGNATURE-----
>

---

From: Anders Logg <[email-masked]>
Date: Fri, 25 Sep 2009 08:22:41 +0200

>   dphi/dz = -constant *phi at zmax

When you integrate by parts you get the term

  -v*dphi/dn*ds

Just replace this by

  -v*dphi/dn*ds = v*(-vphi/dn)*ds = v*constant*phi*ds

This term contains one test function v and one trial function phi, so
just add that term to you bilinear form (it will be in the left-hand
side).

Just stick in a function for the constant that is zero except on z = zmax.

>     void Function::eval(double* values, const double* x, const ufc::cell& ufc_cell,uint cell_index) const;
>
>

This is a low level function call you don't need to worry about.

--
Anders

> i.e. I need to integrate with the trial function as well as the test function
> for this free surface BC?
>
>
> Anyway I have just found that the users used to post to the dofin-dev list and
> there are some similar question to mine there which may be helpful.
> I'll have a dig around in the archives of the dolfin-dev and see if I can work
> it out,
>
>
> regards
>
> Tom
>
>

---

From: Thomas Ward <[email-masked]>
Date: Fri, 2 Oct 2009 18:01:23 +0100

Anders thanks again for your help it is appreciated,

I have tried to put the free surface boundary condition into the bilinear
form,

fs_bc= Function(element)

#fs_bc.eval returns omega^2/9.81 if z=1
a = dot(grad(v), grad(u))*dx +u*v*fs_bc*ds

#g.eval is for the Neuman condtions on x=0 and x=1
L = v*f*dx+v*g*ds

I think I have got that right ?

but when I write the C++ code I have done something wrong because it is not
giving the correct solution (it should generate a free surface gravity wave
but the solution I am getting is as if there is no free surface).
Any ideas?

class Neumann_BC : public Function
{
        void eval ( double* values, const double* x ) const
        {

            if ( x[0] < DOLFIN_EPS ) values[0]= omega;
            else
                if ( x[0]>1-DOLFIN_EPS ) values[0]= -omega;
                else
                     values[0]= 0;
        };
};

class FS_BC : public Function
{
        void eval ( double* values, const double* x ) const
        {
            if ( x[2]>1-DOLFIN_EPS ) values[0]= omega*omega/9.81;
            else values[0]= 0 ;
        };
};

int main ( int argv,char** argc )
{
    int n=10;
    UnitCube mesh ( n, n,n );
    tlaplaceFunctionSpace V ( mesh );
    tlaplaceBilinearForm a ( V, V );
    tlaplaceLinearForm L ( V );

    FS_BC fs_bc;
    a.fs_bc=fs_bc;

    Constant f ( 0.0 );
    L.f = f;

    Neumann_BC g;
    L.g=g;

    std::vector<BoundaryCondition*> bcs;//    bcs.push_back ( &bc0 );
bcs.push_back ( &bc1 );

    VariationalProblem pde ( a, L, bcs );

    Function u;
    pde.solve ( u );

    File file ( "laplace.pvd" );
    file << u;

    return 0;
}

On Fri, Sep 25, 2009 at 7:22 AM, Anders Logg <[email-masked]> wrote:

> work
> >     > this is my first post to this list , and I have just started to try
> to
> >     use ffc
> >     > and dolfin.
> >     > I think I have grasped that to define xmin and xmax neuman
> conditions I
> >     need
> >     > something like
> >     >
> >     > boundary conditions will be enforced naturally by the variational
> form
> >     and it
> >     > I can't see how to imlpenent the zmax boundary condition in dolfin,
> I
> >     think the
> >     >   weak form should be
> >     > no then are there any plans to do so and how difficult would it be
> to do
> >     so?
> >     >
> >     >
> >     > thanks for any help
> >     >
> >     > Tom Ward.
> >
> iEYEARECAAYFAkq8YbEACgkQTuwUCDsYZdGx4gCcD3SRQ+xEjDSWdRiyUqf9xGji
> DjsAn3XMzBUL6oLQjw3e2NDb3VMIEfvy
> =60Wo

---

From: Anders Logg <[email-masked]>
Date: Mon, 5 Oct 2009 19:45:02 +0200

Sorry but I don't have time to debug your model right now... :-)

Let me know if there's anything specific I can answer.

--
Anders

On Fri, Oct 02, 2009 at 06:01:23PM +0100, Thomas Ward wrote:
> I have tried to put the free surface boundary condition into the bilinear form,
>
> giving the correct solution (it should generate a free surface gravity wave but
> the solution I am getting is as if there is no free surface).
> Any ideas?
>

---

From: Thomas Ward <[email-masked]>
Date: Tue, 6 Oct 2009 05:24:50 +0100

No worries Anders, you have helped alot already.

One thing which would help me is some pointers to code where part
boundary integral us carried out in the bilinear form.

  I am reading through some of the papers (when I get time) to get a
better knowledge of the fe method and of the mathematical language
used to describe it and how to pose the physics.   I am also trying to
read through the dolfin source and find sample code of similar
problems. Its a slow process.

Regards

Tom

On 10/5/09, Anders Logg <[email-masked]> wrote:
>>     >     void Function::eval(double* values, const double* x, const
>> ufc::cell&
>>     ufc_cell,uint cell_index) const;
>>     >
>>     >
>>     >   integral_over_surface_of_zmax(constant * trial_function
>> *test_function
>>     )
>>     > ?
>>
>>     > Anyway I have just found that the users used to post to the
>> dofin-dev
>>     list and
>>     >     You will get an integral involving the unknown u on the
>> boundary.
>>     Just
>>     >     > I am impressed with the concept and implementation of your
>> project,
>>     >     thanks you
>>     >     > for making it available.
>>     >     >
>>     >     >
>>     >     > 1.
>>     >     > I am attempting to solve a simple incompressible potential
>> flow
>>     problem
>>     >     > I am slightly worried by this however because I thought that
>> the
>>     dphi/dn
>>     >     > On a separate issue are there any infinite elements
>> implemented in
>>     >     fenics, if

---

From: Anders Logg <[email-masked]>
Date: Tue, 6 Oct 2009 10:13:51 +0200

> boundary integral us carried out in the bilinear form.

# find demo -name '*.ufl' | xargs grep ds
demo/pde/sym-dirichlet-bc/cpp/Poisson.ufl:L = v*f*dx + v*g*ds
demo/pde/lift-drag/cpp/Drag.ufl:M = -p*n[0]*ds
demo/pde/lift-drag/cpp/Lift.ufl:M = p*n[1]*ds
demo/pde/poisson/cpp/Poisson.ufl:L = v*f*dx - v*g*ds
demo/pde/poisson1D/cpp/Poisson.ufl:L = v*f*dx + v*g*ds
demo/pde/dg/advection-diffusion/cpp/AdvectionDiffusion.ufl:a_vel =
dot(jump(v), un('+')*phi('+') - un('-')*phi('-') )*dS  + dot(v,
un*phi)*ds
demo/pde/dg/poisson/cpp/Poisson.ufl:   - dot(grad(v), u*n)*ds \
demo/pde/dg/poisson/cpp/Poisson.ufl:   - dot(v*n, grad(u))*ds \
demo/pde/dg/poisson/cpp/Poisson.ufl:   + (gamma/h)*v*u*ds
demo/pde/elastodynamics/cpp/ElastoDynamics.ufl:# Fields from previous
time step
demo/pde/elastodynamics/cpp/ElastoDynamics.ufl:p       =
Function(element)   # externally applied loads
demo/pde/elastodynamics/cpp/ElastoDynamics.ufl:   + inner(r, f)*dx +
(1.0-alpha_f)*inner(r, p)*ds(3) + alpha_f*inner(r, p0)*ds(3)

> problems. Its a slow process.

You can also look at the preliminary version of the FEniCS book which
has a pretty good introduction FEM. ;-)

http://www.fenics.org/pub/documents/book/fenics-book-2009-06-09-2.pdf

--
Anders

---

From: Thomas Ward <[email-masked]>
Date: Tue, 6 Oct 2009 09:45:22 +0100

Thanks Anders, great book, l particularly like the contents page, it fills
in the gaps between Zienciwicz and C++ code.
Looks like you have already been busy and are going to be even busier!

On Tue, Oct 6, 2009 at 9:13 AM, Anders Logg <[email-masked]> wrote:

> iEYEARECAAYFAkrK/D4ACgkQTuwUCDsYZdFI'�墥짷櫈섵གྷ�瞃
> jNQAmgPcbA3IIJff3ySpoZCoK3wuu6NU
> =apIG