[deal.II] problem with CylindricalManifold<3> after removing cells from a triangulation

[Simon]

Dear all,
 
I created a mesh using GridGenerator::hyper_cube_with_cylindrical_hole and 
removed after that half of the total number of the cells in order to make 
use of symmetry. 
Then I set a SphericalManifold<2> (2D) respectively CylindricalManifold<3> 
(3D) in order to get a true circle (2D) respectively a true cylinder (3D) 
after global and/or adaptive refining. 

Without removing some cells, i.e. just using the 
GridGenerator::hyper_cube_with_cylindrical_hole without using symmetry, 
both 2D and 3D variants work fine. 
With removing some cells, i.e. using symmetry, only in 2D I get a true 
circle as result. In 3D however, the SphericalManifold<3> does not work 
anymore. (I can adaptively refine the inner hole, but the result is not a 
true cylinder anymore which was the case without removing the cells).

In the following is a snippet of my code. It´s hard for me to find out what 
the problem actually is, since in 2D the manifold description worked after 
removing some cells. 

Thanks for helping!

Best
Simon




const double outer_radius = 1.0;
const double inner_radius = 0.5;
const Point center;

GridGenerator::hyper_cube_with_cylindrical_hole(triangulation_tmp,
inner_radius,
 outer_radius,
0.5,
1,
 false 
/*boundary_id_inner_hole is set to 1*/);

std::set::active_cell_iterator> cells_to_remove;
typename Triangulation::active_cell_iterator cell 
=triangulation_tmp.begin_active(),
endc = triangulation_tmp.end();

for(; cell!=endc; cell++)
 {
  if(cell->center()[1]<0)
   {
cells_to_remove.insert(cell);
}
  }
GridGenerator::create_triangulation_with_removed_cells(triangulation_tmp, 
cells_to_remove, triangulation);

std::unique_ptr> ptr_manifold=nullptr;

if(dim==2)
{
ptr_manifold = std::make_unique>(center);
}
else if(dim==3)
{
ptr_manifold = std::make_unique>(dim-1);
}
else
{
throw std::runtime_error("only allowed for dim == 2 or dim == 3");
}
types::boundary_id boundary_id_inner_hole=1;
types::manifold_id manifold_id_inner_hole=1;
   
triangulation.set_all_manifold_ids_on_boundary(boundary_id_inner_hole,manifold_id_inner_hole);

triangulation.set_manifold (manifold_id_inner_hole, *ptr_manifold);

..

-- 
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/ab270946-bbe3-4952-a862-b7fd2e6ccd3bn%40googlegroups.com.

Re: [deal.II] interpolate FE_Nelelec

[Konrad Simon]

Hi John,

Maybe I can give some hints of how I would approach this problem (these are 
just some quick thoughts):
Write this problem as a vector Laplace problem

curl(nu*curl(A)) - grad(div(A)) = J

which you then have to write as a system of PDEs. Note, this can only be 
the same as

curl(nu*curl(A)) = J

if J is a curl itself.

If you do not want to impose conditions on the divergence then you have two 
options:

1.) You use sigma = div(A) as an auxiliary variable and impose natural 
boundary conditions (i.e., the ones that you get through integration by 
parts) on A*n and nu*curl(A)xn, where n is the outward unit normal. This 
means that sigma is a 0-form and A is a 1-form -> sigma is then discretized 
(for example) by FE_Q(n+1) and A by FE_Nedelec(n)

System:

(A1) sigma + div(A) = 0
(B1) grad(sigma) + curl(nu*curl(A)) = J

For the weak form you multiply (A1) with a FE_Q test function and integrate 
the second summand by parts. Secondly, you multiply (B1) with a Nedelec 
test function and also integrate the second summand by parts. You will see 
your natural boundary conditions pop up.

2.) You use sigma = nu*curl(A) as an auxiliary variable and impose 
essential boundary conditions on A*n and nu*curl(A)xn, where n is the 
outward unit normal. This means that sigma is interpreted as a 1-form and A 
is a 2-form -> sigma is then discretized (for example) by FE_Nedelec(n) and 
A by FE_RaviartThomas(n)

System:

(A2) (1/nu)*sigma - curl(A) = 0
(B2) curl(sigma) -  grad(div(A)) = J

For the weak form you multiply (A2) with a Nedelec test function and 
integrate the second summand by parts. Secondly, you multiply (B2)
with a Raviart-Thomas test function and also integrate the second summand 
by parts. You will NOT see your natural boundary conditions pop up since 
conditions on A*n and nu*curl(A)xn are essential in this case. You need to 
enforce them on the function spaces directly. In deal.ii you do so by using 
project_boundary_values_curl_conforming_l2 

 
or project_boundary_values_div_conforming 


Note that the additional boundary conditions make the system invertible and 
if J is a curl it will turn out that div(A)=0.

There is a field in mathematics that is called finite element exterior 
calculus (FEEC) that answers the question of stability when solving such 
problems. See this book 
. I guess Chapters 
4, 5 and 6 are most interesting for you.

> I do not understand when a geometry gets complicated. Is a toroid inside a 
> sphere, both centred at the origin, a complicated geometry? To start with, 
> I will use a relatively simple mesh. I will not use local refinement, only 
> global. I can do without refinement as well, i.e. making the mesh in gmsh.
>
> Yes, I would like to know more about orientations issues on complicated 
> geometries. Could please tell me about it? I would very much prefer to use 
> FE_Nedelec without hierarchical shape functions. The first order 3D 
> FE_Nedelec will do nicely provided the orientation issues will be 
> irrelevant to the mesh.
>
 Lowest order Nedelec and all other elements I mentioned are fine on the 
meshes you need I guess.

Hope that helps.
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/c2773841-4e03-42a1-9f1e-e1a41c0933c9n%40googlegroups.com.

[deal.II] integrating over boundaries

[Sylvain Mathonnière]

Hello everyone !

I have a very basic question but somehow I did not manage to find the 
answer in tutorials, neither looking at the documentation or in the this 
user group. Or more probably I did not understand it when I saw it since I 
am fairly new to deal.II and to C++.

*My main goal :* I am trying to integrate over a boundary (a line since I 
am in 2D) the heat flux to get the total heat flux through my boundary.

*What I did :*
1 - I looked at this documentation 
https://www.dealii.org/current/doxygen/deal.II/classDataPostprocessorVector.html
 
to calculate first the heat flux on my boundary which worked fine.

2 - I tried to integrate the norm of the heat flux over my boundary in the 
function *evaluate_scalar_field()* function but I could not. My best effort 
is this :

template 
class HeatFluxPostprocessor : public DataPostprocessorVector
{
public:
  HeatFluxPostprocessor ()
:
DataPostprocessorVector ("_k_grad_T",
  update_gradients | update_quadrature_points)
  {}
  virtual void evaluate_scalar_field(
const DataPostprocessorInputs::Scalar _data,
std::vector >   _quantities) const
{
  // I recover the value of conductivity
  MaterialData material_data;
  const typename DoFHandler::cell_iterator current_cell = 
input_data.template get_cell>();
  FullMatrix  conductivity = material_data.get_th_conductivity(
current_cell->material_id());

  double heat_flux_cell_norm = 0.;

  for (auto face_index : GeometryInfo::face_indices())
  {
  if (current_cell->face(face_index)->boundary_id() == 2)
  {
double long_cell = 0;
long_cell = (current_cell->face(face_index)->vertex(0) - 
current_cell->face(face_index)->vertex(1)).norm();

for (unsigned int d=0; dhttp://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/6058222b-5470-4bbd-ad85-7f43327b37adn%40googlegroups.com.