Re: [deal.II] Boundary conditions on components in with FESystem

[Konrad Simon]

Dear Wolfgang,

I may have stumbled over something. I created a 3D hyper_shell with 6 
coarse cells which is nothing but a cube with a cubical void in the center. 
I then asked for the face orientations of each face in each cell. In each 
of the 6 cells all faces have standard orientation according to the deal.ii 
convention. Mathematically this is not possible if you want a consistent 
orientation of normals cross the mesh (right?). I think, two cells must 
have flipped two flipped face normals for consistency.

In the GeometryInfo documentation it says:

However, it turns out that a significant number of 3d meshes cannot satisfy 
this convention. This is due to the fact that the face convention for one 
cell already implies something for the neighbor, since they share a common 
face and fixing it for the first cell also fixes the normal vectors of the 
opposite faces of both cells. It is easy to construct cases of loops of 
cells for which this leads to cases where we cannot find orientations for 
all faces that are consistent with this convention.

For this reason, above convention is only what we call the *standard 
orientation*. deal.II actually allows faces in 3d to have either the 
standard direction, or its opposite, in which case the lines that make up a 
cell would have reverted orders, and the above line equivalences would not 
hold any more. You can ask a cell whether a given face has standard 
orientation by calling cell->face_orientation(face_no): if the result 
is true, then the face has standard orientation, otherwise its normal 
vector is pointing the other direction. There are not very many places in 
application programs where you need this information actually, but a few 
places in the library make use of this. Note that in 2d, the result is 
always true. More information on the topic can be found in this glossary 

 article.
Is there a simple way to tell deal.ii to flip the normal of such a face? If 
yes, does this affect the edge orientation? I am asking since this is 
important for mapping Nedelec shape functions and I do not want to solve 
the problem for faces create one for edges.

A way out could be to simply flip the signs of (mapped) Raviart-Thomas 
shape functions that are supported on a non-consistent face. I guess this 
strategy would work for any H(div)-conformal shape function. Am I correct?

I am willing to fix that (in deal.ii itself) if necessary. 

Best,
Konrad

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/fc6bebbb-7f36-4e63-9028-f3ced9ba82f8n%40googlegroups.com.

Re: [deal.II] Boundary conditions on components in with FESystem

[Wolfgang Bangerth]

On 11/6/20 10:11 AM, Konrad Simon wrote:


I have a code version in which I am using a FESystem (3d) composed of three 
elements: FE_Nedelec (w), FE_RaviartThomas (u) and DGQ (p).


The code compiles fine but gives unreasonable results and I do not see what I 
am doing wrong or what I forgot.
Instead of trying to divine what the code is doing (and apparently doing 
wrong), let me outline the way I would think about this:


* You say that the results are "unreasonable". That can mean a lot of 
different things. Try to explain for yourself what you see, and understand 
what part looks wrong. Is it the boundary values? That is, does the solution 
you see visualized violate the boundary conditions you wanted to impose? If 
so, that's clearly an area of concern.


* If the solution does seem to satisfy the boundary conditions, but is wrong 
in the interior of the domain, then the problem lies elsewhere.


* In any case, if you can't find where things go wrong, try solving a 
decoupled problem where you assemble a matrix in which each variable appears 
in only one equation. This helps you identify *which boundary condition*, or 
*which equation* is wrongly implemented.


In the end, debugging comes down to generating hypotheses about what could be 
wrong, and then testing each hypothesis in detail until you can confirm or 
deny that your hypothesis is the cause of the problem. Creating hypotheses 
often comes down to carefully looking at the solution and verbalizing for 
yourself what you see and how it (qualitatively) differs from what you expect. 
In practice, this often leads to faster ways to debugging than randomly poking 
at individual pieces of your code.


Best
 W.


--

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

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/e41c1a04-c857-ab3f-e2d8-6c9178fdaeb6%40colostate.edu.

[deal.II] Boundary conditions on components in with FESystem

[Konrad Simon]

Dear all,

I have a code version in which I am using a FESystem (3d) composed of three 
elements: FE_Nedelec (w), FE_RaviartThomas (u) and DGQ (p). 

The code compiles fine but gives unreasonable results and I do not see what 
I am doing wrong or what I forgot.

I need to impose (essential) boundary conditions on the first two 
components: 

   - on w: w x n =0 (n is normal)
   - on u: u*n=0 (no normal flux)

I assume that I am doing something wrong in the setup of BCs or the system 
(I guess the solver is likely ruled since two different solvers give 
exactly the same nonsense). This is the code snippet that I suspect to be 
buggy. I hope the names are self explanatory (essentially it is a 
modification of step-32). Note that dim=3 here.

Any help would be appreciated. 

Best,
Konrad

/
  // System and dof setup
  /

  template 
  void
  BoussinesqModel::setup_nse_matrices(
const std::vector _partitioning,
const std::vector _relevant_partitioning)
  {
nse_matrix.clear();
LA::BlockSparsityPattern sp(nse_partitioning,
nse_partitioning,
nse_relevant_partitioning,
this->mpi_communicator);

Table<2, DoFTools::Coupling> coupling(2 * dim + 1, 2 * dim + 1);
for (unsigned int c = 0; c < 2 * dim + 1; ++c)
  {
for (unsigned int d = 0; d < 2 * dim + 1; ++d)
  {
if (c < dim)
  {
if (d < 2 * dim - 1)
  coupling[c][d] = DoFTools::always;
else
  coupling[c][d] = DoFTools::none;
  }
else if ((c >= dim) && (c < 2 * dim))
  {
coupling[c][d] = DoFTools::always;
  }
else if (c == 2 * dim)
  {
if ((d >= dim) && (d < 2 * dim))
  coupling[c][d] = DoFTools::always;
else
  coupling[c][d] = DoFTools::none;
  }
  }
  }

DoFTools::make_sparsity_pattern(nse_dof_handler,
coupling,
sp,
nse_constraints,
false,
Utilities::MPI::this_mpi_process(
  this->mpi_communicator));
sp.compress();

nse_matrix.reinit(sp);
  }

  template 
  void
  BoussinesqModel::setup_dofs()
  {
TimerOutput::Scope timing_section(
  this->computing_timer, "BoussinesqModel - setup dofs of systems");

/*
 * Setup dof handlers for nse and temperature
 */
std::vector nse_sub_blocks(3, 0);
nse_sub_blocks[1] = 1;
nse_sub_blocks[2] = 2;

nse_dof_handler.distribute_dofs(nse_fe);

/*
 * Reduce bandwidth for ILU preconditioning
 */
DoFRenumbering::Cuthill_McKee(nse_dof_handler);

DoFRenumbering::block_wise(nse_dof_handler);

/*
 * Count dofs
 */
std::vector nse_dofs_per_block(3);
DoFTools::count_dofs_per_block(nse_dof_handler,
   nse_dofs_per_block,
   nse_sub_blocks);
const unsigned int n_w = nse_dofs_per_block[0], n_u = 
nse_dofs_per_block[1],
   n_p = nse_dofs_per_block[2],
   n_T = temperature_dof_handler.n_dofs();

/*
 * Comma separated large numbers
 */
std::locale s = this->pcout.get_stream().getloc();
this->pcout.get_stream().imbue(std::locale(""));

/*
 * Print some mesh and dof info
 */
this->pcout << "Number of active cells: "
<< this->triangulation.n_global_active_cells() << " (on "
<< this->triangulation.n_levels() << " levels)" << std::endl
<< "Number of degrees of freedom: " << n_w + n_u + n_p + n_T
<< " (" << n_w << " + " << n_u << " + " << n_p << " + " << 
n_T
<< ")" << std::endl
<< std::endl;
this->pcout.get_stream().imbue(s);

/*
 * Setup partitioners to store what dofs and matrix entries are stored 
on
 * the local processor
 */
{
  nse_index_set = nse_dof_handler.locally_owned_dofs();

  nse_partitioning.push_back(nse_index_set.get_view(0, n_w));
  nse_partitioning.push_back(nse_index_set.get_view(n_w, n_w + n_u));
  nse_partitioning.push_back(
nse_index_set.get_view(n_w + n_u, n_w + n_u + n_p));

  DoFTools::extract_locally_relevant_dofs(nse_dof_handler,
  nse_relevant_set);

  nse_relevant_partitioning.push_back(nse_relevant_set.get_view(0, 
n_w));
  nse_relevant_partitioning.push_back(
nse_relevant_set.get_view(n_w, n_w + n_u));
  nse_relevant_partitioning.push_back(