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 <[EMAIL PROTECTED]>:
>>>>>
>>>>>
>>>>>> 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.
/murtazo
_______________________________________________
DOLFIN-dev mailing list
[email protected]
http://www.fenics.org/mailman/listinfo/dolfin-dev
