Representing a void subdomain with a refined triangulation

[Alexander]

Dear all,

I have a triangulation which has been refined along some (polygonal) patch, thus forming two non-overlapping domains (say, domain #1 and #2). My task is to extract cells for domain #1, and solve BVP for this domain (with known boundary conditions). The rest of the cells (domain #2) are not needed. An example of such situation is shown attached (with land and ocean/sea being two domains).

FWIW, the problem is posed in dim = 2, spacedim = 3. The solution is in H1.

Here is what I struggle with:

- I thought of using GridGenerator::create_triangulation_with_removed_cells, but this one does not accept refined non-conforming meshes.

- My next thought was to detect those edges, which domains #1 and #2 share and calculate face integrals on them during assembly, thus incorporating BCs. Will DoFTools::make_flux_sparsity_pattern be the right function for this?

- Next, what should I do on cells in domain #2? Logically, I should somehow treat them with FE_Nothing... But how to get this done right?

Thanks!

Alexander

[Wolfgang Bangerth]

Alexander,

> I have a triangulation which has been refined along some (polygonal) patch,
> thus forming two non-overlapping domains (say, domain #1 and #2). My task is
> to extract cells for domain #1, and solve BVP for this domain (with known
> boundary conditions). The rest of the cells (domain #2) are not needed. An
> example of such situation is shown attached (with land and ocean/sea being two
> domains).
>
> FWIW, the problem is posed in dim = 2, spacedim = 3. The solution is in H1.
>
> Here is what I struggle with:
>
> - I thought of using GridGenerator::create_triangulation_with_removed_cells,
> but this one does not accept refined non-conforming meshes.

Yes. You can't remove cells this way because then you would have (non-active)
cells for which not all children are part of the mesh. That is not allowed.

> - My next thought was to detect those edges, which domains #1 and #2 share and
> calculate face integrals on them during assembly, thus incorporating BCs. Will
> DoFTools::make_flux_sparsity_pattern be the right function for this?

You mean you want to impose boundary conditions weakly? Like through the
Nitsche method? You could do that. But you can also just impose them strongly.
In essence, duplicate the functionality of the
VectorTools::interpolate_boundary_values() function but instead of running it
on only the faces at the boundary of the mesh, you run it also on the faces at
the boundaries between the two subdomains.

> - Next, what should I do on cells in domain #2? Logically, I should somehow
> treat them with FE_Nothing... But how to get this done right?

If you impose Dirichlet boundary conditions on the interface between the
subdomains, then no knowledge of domain #2 leaks into domain #1. In other
words, you can do whatever you please on domain #2. One option would be to
just solve the same equation there -- you'd get something for the solution
there, but you really don't care very much what exactly it is :-) If you
wanted to reduce the numerical effort, you could just assemble a mass matrix
on domain #2 instead of a stiffness matrix.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Alexander]

Thanks, Wolfgang!

What you describe makes perfect sense in case of Dirichlet BCs. However, I have no-flux (natural) BCs at the boundary between domains. I don't understand what to do with the cells in domain #2 to satisfy these?

Additionally, would it be possible to further reduce the solution effort by constraining dofs inside domain #2 to zero?

Best

Alexander

[Wolfgang Bangerth]

> What you describe makes perfect sense in case of Dirichlet BCs. However, I
> have no-flux (natural) BCs at the boundary between domains. I don't understand
> what to do with the cells in domain #2 to satisfy these?

Hm, I must admit that I don't quite know what the right solution is then. It's
*possible* that just adding the Neumann boundary terms on the interfaces is
enough. But I'm not entirely sure and it would seem to me that what happens in
domain #2 then still influences what is going on in domain #1 -- the Neumann
terms are just another source term, and the ellipticity of the problem in that
case implies that whatever sources you have in domain #2 (or don't have) have
an effect on the solution in domain #1. You might want to think this through
some more.

This would suggest that using FE_Nothing in domain #2 might be the better
choice after all. step-46 will be your friend understanding how to use it.

I'll add another thought: If you have non-zero Neumann boundary conditions,
and you approximate the original domain by one that approximates the boundary
by a stair-step function as you seem to do, then the numerical solution does
not converge to the exact one. That's because the stair-step boundary has a
different length than the exact boundary (and doesn't converge to the exact
length), and the total flux into the domain equals the average flux on the
boundary times the length of the boundary.

(Fun fact: Ivo Babuska explained this to me at a conference dinner at the
first conference I ever went to, in 1998 or 1999.)

> Additionally, would it be possible to further reduce the solution effort by
> constraining dofs inside domain #2 to zero?

Yes, but that's equivalent to using a diagonal mass matrix. I suspect it's not
worth it.

[Alexander]

Wolfgang,

thanks, I'm inclined towards using FE_Nothing. Since this becomes quite technical, I'm thinking of a following test. Can you please confirm if this is a sensible way of validation.

(i) Construct a single cell quad mesh. Possibly refine n times.

(ii) At the top (arbitrary choice) boundary set no-flux BCs. Set (smooth) Dirichlet BCs on all other sides.

(iii) Solve laplace problem.

(iv) Construct a two cells quad mesh such that the top boundary of the mesh from (i) became an internal face shared by the two cells.

(v) Set cell #1 to FE_Q and cell #2 to FE_Nothing. Refine n times.

(vi) Impose no-flux face integral at the shared faces. Set same Dirichlet BCs as in (ii) on all other sides of cell #1.

(vii) Solve laplace problem.

(viii) Compare solution in cell #1 to the one in step (iii). They should be equal.

I'll add another thought: If you have non-zero Neumann boundary conditions,
and you approximate the original domain by one that approximates the boundary
by a stair-step function as you seem to do, then the numerical solution does
not converge to the exact one. That's because the stair-step boundary has a
different length than the exact boundary (and doesn't converge to the exact
length), and the total flux into the domain equals the average flux on the
boundary times the length of the boundary.

An interesting point, although I guess in a row of all approximations I have already done in this problem, this one won't be the first to hit me back.

[Wolfgang Bangerth]

On 6/18/19 2:35 PM, Alexander wrote:
> Wolfgang,
> thanks, I'm inclined towards using FE_Nothing. Since this becomes quite
> technical, I'm thinking of a following test. Can you please confirm if this is
> a sensible way of validation.
>
> (i) Construct a single cell quad mesh. Possibly refine n times.
> (ii) At the top (arbitrary choice) boundary set no-flux BCs. Set (smooth)
> Dirichlet BCs on all other sides.
> (iii) Solve laplace problem.
>
> (iv) Construct a two cells quad mesh such that the top boundary of the mesh
> from (i) became an internal face shared by the two cells.
> (v) Set cell #1 to FE_Q and cell #2 to FE_Nothing. Refine n times.
> (vi) Impose no-flux face integral at the shared faces. Set same Dirichlet BCs
> as in (ii) on all other sides of cell #1.
> (vii) Solve laplace problem.
>
> (viii) Compare solution in cell #1 to the one in step (iii). They should be equal.

Yes, this seems eminently sensible.

> I'll add another thought: If you have non-zero Neumann boundary conditions,
> and you approximate the original domain by one that approximates the boundary
> by a stair-step function as you seem to do, then the numerical solution does
> not converge to the exact one. That's because the stair-step boundary has a
> different length than the exact boundary (and doesn't converge to the exact
> length), and the total flux into the domain equals the average flux on the
> boundary times the length of the boundary.
>
>
> An interesting point, although I guess in a row of all approximations I have
> already done in this problem, this one won't be the first to hit me back.

It should only affect you if you had non-zero boundary fluxes. But it's an
O(1) error, and it's not small because it is proportional to the ratio of the
lengths of the exact boundary to the length of the approximated stair-step
boundary. If your exact domain is a triangle, then that ratio is sqrt(2)=1.4!

[Alexander]

It should only affect you if you had non-zero boundary fluxes. But it's an
O(1) error, and it's not small because it is proportional to the ratio of the
lengths of the exact boundary to the length of the approximated stair-step
boundary. If your exact domain is a triangle, then that ratio is sqrt(2)=1.4!

Oh, I see... Fortunately, I have zero flux conditions at the sea/land interface.

The domain boundaries have inhomogeneous Dirichlet BCs, but they are straight anyway.

Thanks for your help. I will try this simple test and see how it goes.

Best

Alexander