dolfin team mailing list archive

["Murtazo Nazarov" <murtazo@xxxxxxxxxxx>]

> 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

>
> --
> Anders
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/dolfin-dev
>