dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Tue, Apr 01, 2008 at 11:37:58PM +0200, Murtazo Nazarov wrote:
> > On Tue, Apr 01, 2008 at 09:52:02PM +0200, Murtazo Nazarov wrote:
> >> > Quoting Murtazo Nazarov <murtazo@xxxxxxxxxxx>:
> >> >
> >> >> Anders Logg wrote:
> >> >> > On Tue, Apr 01, 2008 at 02:32:17PM +0200, Murtazo Nazarov wrote:
> >> >> >
> >> >> >> Anders Logg wrote:
> >> >> >>
> >> >> >>> On Mon, Mar 31, 2008 at 07:15:23PM +0200, Murtazo Nazarov wrote:
> >> >> >>>
> >> >> >>>
> >> >> >>>> Anders Logg wrote:
> >> >> >>>>
> >> >> >>>>
> >> >> >>>>> On Mon, Mar 31, 2008 at 03:52:57PM +0200, Murtazo Nazarov
> >> wrote:
> >> >> >>>>>
> >> >> >>>>>
> >> >> >>>>>
> >> >> >>>>>> Hi,
> >> >> >>>>>>
> >> >> >>>>>> I want to apply a boundary condition to the vertices on a
> >> >> boundary.
> >> >> For that I define:
> >> >> >>>>>>
> >> >> >>>>>> // Sub domain for MyBC
> >> >>
> >> >> >>>>>> class MyBC_Boundary2D : public SubDomain
> >> >> >>>>>> {
> >> >> >>>>>> public:
> >> >> >>>>>>   bool inside(const real* p, bool on_boundary) const
> >> >> >>>>>>   {
> >> >> >>>>>>     return  on_boundary && (p[0] < xmax - bmarg) && (p[0] >
> >> xmin
> >> >> +
> >> >> bmarg);
> >> >> >>>>>>   }
> >> >> >>>>>> };
> >> >> >>>>>>
> >> >> >>>>>> Then I initialize:
> >> >> >>>>>>
> >> >> //-----------------------------------------------------------------------------
> >> >>
> >> >> >>>>>> void MyBC::init(SubDomain& sub_domain)
> >> >> >>>>>> { ...
> >> >> >>>>>>   mesh.init(0);
> >> >> >>>>>>   sub_domains = new MeshFunction<uint>(mesh, 0);
> >> >> >>>>>> ...}
> >> >> >>>>>>
> >> >> >>>>>>
> >> >> >>>>>>
> >> >> >>>>> Do you remember to set everything to 1 (number of subdomains)
> >> >> here?
> >> >> >>>>>
> >> >> >>>>>
> >> >> >>>>>
> >> >> >>>>>
> >> >> >>>> Yes I do:
> >> >> >>>>
> >> >> >>>>   // Mark everything as sub domain
> >> >> >>>> 1
> >> >> >>>>   (*sub_domains) = 1;
> >> >> >>>>
> >> >> >>>>
> >> >> >>>>   // Mark the sub domain as sub domain
> >> >> >>>> 0
> >> >> >>>>   sub_domain.mark(*sub_domains, 0);
> >> >> >>>>
> >> >> >>>> /murtazo
> >> >> >>>>
> >> >> >>>>
> >> >> >>> Then I have no idea. You just need to dig into the code and see
> >> what
> >> >> >>> goes wrong.
> >> >> >>>
> >> >> >>>
> >> >> >>>
> >> >> >> I think I found a mistake in Subdomain.mark(). The thing is bool
> >> >> >> on_boundary should be defined
> >> >> >>  inside the loop for computing subdomain marker:
> >> >> >>
> >> >> >>   // Always false when not marking
> >> >> >> facets
> >> >>
> >> >> >>
> >> >> >>   //bool on_boundary =
> >> >> >> false;
> >> >>
> >> >> >>
> >> >> >>
> >> >> >>   // Compute sub domain
> >> >> >> markers
> >> >>
> >> >> >>
> >> >> >>   for (MeshEntityIterator entity(mesh, dim); !entity.end();
> >> ++entity)
> >> >> >>   {
> >> >> >>     bool on_boundary = false;
> >> >> >>     ...
> >> >> >>    }
> >> >> >>
> >> >> >> Otherwise, if you want to apply a boundary condition on vertices,
> >> it
> >> >> >> becomes always
> >> >> >> true after once meeting a vertex on the boundary.
> >> >> >>
> >> >> >>
> >> >> >> It solved problem for me.
> >> >> >>
> >> >> >
> >> >> > I'm not sure. If you are marking vertices, then they will never be
> >> on
> >> >> > the boundary, unless you are doing this for a 1-dimensional mesh
> >> >> > (which typically has just 2 boundary vertices). Only *facets* will
> >> >> > ever be marked as being on the boundary, so for a 2-dimensional
> >> mesh,
> >> >> > edges on the boundary will be marked as being on the boundary and
> >> for
> >> >> > a 3-dimensional mesh, faces on the boundary will be marked as being
> >> on
> >> >> > the boundary.
> >> >> >
> >> >> >
> >> >> With the change I wrote I was able to mark vertices. Because,
> >> >> if I mark facets or faces, I miss some nodes in the boundary
> >> >> of the subdomain where I am applying my bc. That destroys
> >> >> the solution. Now, I apply my boundary only for vertices
> >> >> I need, and I do not miss any nodes. It work perfectly with
> >> >> the slip BC.
> >> >
> >> > Sounds like a similar issue that we have for Discontinuous elements.
> >> Have
> >> > you
> >> > tried to apply BCs using a geometrical (or pointwise) approach?
> >> >
> >>
> >> Yes, I have. But as I said I had a problem with some side nodes. The
> >> boundary condition did not apply to those nodes. Therefore I got
> >> incorrect
> >> result. But, with this approach everything went well.
> >>
> >> /murtazo
> >
> > I still don't understand what the problem is. What kind of mesh do you
> > have? Triangles or tetrahedra? What are you marking? Vertices?
> >
> > If you are marking vertices, then on_boundary should always be
> > false. It is set to false on line 53 and it will never change since
> > line 60 will not be reached.
> 
> Ok, I see. It seems we are looking to some different version of DOLFIN. My
> Subdomain.init looks like the following:
> 
>   // Always false when not marking facets
>   //bool on_boundary = false;
> 
>   // Compute sub domain markers
>   for (MeshEntityIterator entity(mesh, dim); !entity.end(); ++entity)
>   {
>     // on_boundary should be initialized here:
>     bool on_boundary = false;
> 
>     // Check if entity is on the boundary if entity is a facet
>     if (dim == D - 1)
>       on_boundary = entity->numEntities(D) == 1;
>     if (dim == 0)
>     {
>       for (FacetIterator fi(*entity); !fi.end(); ++fi)
>       {
>         if(fi->numEntities(D) == 1)
>           on_boundary = true;
>       }
>     }
>     ...
>    }
> 
> Now, you see then in the loop FacetIterator on_boundary is true, even if
> we iterate over the vertices. Here dim = 0.
> 
> /murtazo

ok, good. So does it work for you with the latest version?

-- 
Anders