dolfin team mailing list archive

[Anne Voigt <anne.voigt@xxxxxxxxxxxxxxxxxxxxx>]

Am Dienstag, den 10.06.2008, 19:06 +0200 schrieb Anders Logg:
> On Fri, Jun 06, 2008 at 10:28:54AM +0200, Anne Voigt wrote:
> > Am Dienstag, den 03.06.2008, 22:02 +0200 schrieb Anders Logg:
> > > On Fri, May 30, 2008 at 10:02:29AM +0200, Anne Voigt wrote:
> > > > Hello all,
> > > > 
> > > > I'm working with FFC for quite a while and I like the easy way to
> > handle
> > > > multilinear forms with the help of FFC.
> > > > 
> > > > But now I have a problem where I'm not sure if it is possible to
> > > > programme this with FFC.
> > > > 
> > > > I'm working on an optimal control problem where I have to solve an
> > > > integral over just a part of the whole twodimensional domain (x,y)
> > \in
> > > > \Omega. So what I need is the bilinear form
> > > > 
> > > > \int_{y_0}^{y_1} u*v dy where [y_0,y_1] \in \Omega
> > > > 
> > > > Is it possible to implement something like that in FFC and how can I
> > do
> > > > that?  
> > > > 
> > > > Hope somebody can help me?!
> > > 
> > > Yes, that's easy (but not so well documented I admit).
> > > 
> > > The simplest way is to just supply a SubDomain as argument to the
> > assembler:
> > > 
> > >   assemble(A, form, mesh, sub_domain);
> > > 
> > > Create a SubDomain which overloads the inside() function, returning
> > > True when you are inside your subdomain, for example
> > > 
> > >   class MySubDomain : public SubDomain
> > >   {
> > >     bool inside(const real* x, bool on_boundary) const
> > >     {
> > >       return x[0] + x[1] > 1.0;
> > >     }
> > >   };  
> > > 
> > > If you can't specify your domain based on an expression, you need to
> > > create a MeshFunction that labels which cells are included in your
> > > domain.
> > > 
> > > It's also possible to have different integrals over different
> > > subdomains, integrals over subsets of the boundary etc.
> > > 
> > 
> > Thanks for your respond. But it seems like it would not work.
> > 
> > If I do
> > 
> > #include <iostream>
> > #include <fstream>
> > #include <sstream>
> > #include <vector>
> > 
> > #include <dolfin.h>
> > #include <dolfin/dolfin_mf.h>
> > 
> > 
> > #include "MassMatrix2D.h"
> > 
> > using namespace dolfin;
> > using namespace std;
> > 
> > 
> > class MySubDomain : public SubDomain
> > {
> >         bool inside(const dolfin::real *x, bool on_boundary) const
> >         {       
> >         return true;    
> >         }
> > };
> > 
> > int main(int argc, char* argv[])
> > {
> >   dolfin_init(argc, argv);
> >   Mesh mesh("cylinder_10.xml");
> > 
> >   int inode = mesh.numVertices();
> >   PETScMatrix K;
> >   MassMatrix2DBilinearForm a;
> >   MySubDomain subDomain;
> >   assemble(K, a, mesh, subDomain);
> > }
> > 
> > I get there error message
> > 
> > | Assembling over cells                                           |
> > |-----------------------------------------------------------------| 0.0%
> > [0]PETSC ERROR:
> > ------------------------------------------------------------------------
> > [0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation,
> > probably memory access out of range
> > [0]PETSC ERROR: Try option -start_in_debugger or
> > -on_error_attach_debugger
> > [0]PETSC ERROR: or see
> > http://www.mcs.anl.gov/petsc/petsc-as/documentation/troubleshooting.html#Signal[0]PETSC ERROR: or try http://valgrind.org on linux or man libgmalloc on Apple to find memory corruption errors
> > [0]PETSC ERROR: likely location of problem given in stack below
> > [0]PETSC ERROR: ---------------------  Stack Frames
> > ------------------------------------
> > [0]PETSC ERROR: Note: The EXACT line numbers in the stack are not
> > available,
> > [0]PETSC ERROR:       INSTEAD the line number of the start of the
> > function
> > [0]PETSC ERROR:       is given.
> > [0]PETSC ERROR: --------------------- Error Message
> > ------------------------------------
> > [0]PETSC ERROR: Signal received!
> > [0]PETSC ERROR:
> > ------------------------------------------------------------------------
> > [0]PETSC ERROR: Petsc Release Version 2.3.3, Patch 8, Fri Nov 16
> > 17:03:40 CST 2007 HG revision:
> > 414581156e67e55c761739b0deb119f7590d0f4b      
> > ...
> > 
> > But I should get the same result as I would get if I woulde assemble
> > over the whole domain or???
> > It does also not work if would compile it with your little example. So
> > is it possible that there is a bug in one of the functions?!
> > 
> > Anne                                                    
> 
> Works for me. Take a look in
> 
>   sandbox/misc/
> 
> Does this give you a segmentation fault?

I think there was a misunderstanding. I'm able to assemble the Stiffness
and Mass Matrix over the whole domain. 
But it doesn't work for any(!) subdomain.

I had a look in

sandbox/misc

but there the whole Stiffness Matrix is assembled (and not just a part).
This works fine, but as soon as I use 

 assemble(K, a, mesh, subDomain);

I get a segmentation fault. So first I thought was that my
implementation of the subDomain was wrong. So I tried 

 class MySubDomain : public SubDomain
{         bool inside(const dolfin::real *x, bool on_boundary) const
        {       
         return true;    
         }
 };

Now I thought I would get the WHOLE Stiffness Matrix (or what ever I had
as a ffc form). But I got a segmentation fault again. How is that
possible?

If you like you can try it by yourself. Here is my main-program 


#include <dolfin.h>


#include "MassMatrix2D.h"

using namespace dolfin;

class MySubDomain : public SubDomain
{
	bool inside(const dolfin::real *x, bool on_boundary) const
	{	
	return true;	
	}
};

int main(int argc, char* argv[])
{
  dolfin_init(argc, argv);
  UnitSquare mesh(10, 10);

  int inode = mesh.numVertices();
  PETScMatrix K;
  MassMatrix2DBilinearForm a;
  MySubDomain subDomain;
  assemble(K, a, mesh, subDomain);

  return 0;
}


You just have to include the MassMatrix2D.h ( src/kernel/mf/ffc-forms ).
I think there is something wrong in the assemble function but of course
I'm not sure...

It would be great if you would be able to help!

Thanks Anne