dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Wed, Apr 02, 2008 at 01:01:41PM +0200, Murtazo Nazarov wrote:
> Anders Logg wrote:
> > 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?
> >
> >   
> Yes, It works for me with the latest version of Subdomain.cpp (from 
> dolfin-dev) and with dolfin-snapshot-2008-01-21.

Excellent! Let that be a lesson to you, to always use the latest
version. ;-)

-- 
Anders