Re: [deal.II] Averaging derived solution quantities (mises stress)

[Konrad Schneider]

estion of the FAQ the cell-loop then has to 
start like this:
for(const auto & tria_cell : this->triangulation.active_cell_iterators()){ 
typename DoFHandler::active_cell_iterator cell (>triangulation, 
tria_cell->level(), tria_cell->index(), _handler_recovery);
typename DoFHandler::active_cell_iterator stress_cell (>
triangulation, tria_cell->level(), tria_cell->index(), >dof_handler);
However, in contrast to the FAQ answer one has to use the "typename" 
keyword.

Many thanks again and best regards
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/edbc2f13-176f-47f0-85e5-e3a0bd625ab0n%40googlegroups.com.

Re: [deal.II] Averaging derived solution quantities (mises stress)

[Konrad Schneider]

there was a small bug in my code I posted above. the following code is 
running fine.

Konrad Schneider schrieb am Freitag, 17. März 2023 um 17:03:18 UTC+1:

> Dear Wolfgang,
> yes the mesh is quite coarse, but I want to investigate the functionality 
> of how deal.ii outputs stresses and how these stresses are then displayed 
> via for instance paraview. To really see the effects I intentionally use 
> the coarse mesh.
> So far, I actually do create a (derived) DataPostprocessor-Object (see 
> the thinned out and altered step-8 tutorial attached to this message) 
> similar to the procedure in step-32. You wrote, that the gradients are 
> computed at the vertices. How so? If I look at line 291 of my code or at 
> step-32, there is a loop over all quadrature points (  for (unsigned int q 
> = 0; q < n_quadrature_points; ++q) ) which gets called for every cell. In 
> the scope of this loop the shape function gradients at the quadrature 
> points are used to compute the desired output variable by a special 
> calculation rule. However, these are not the vertices. If the gradients are 
> somehow again computed at the vertices, which might be an internal 
> procedure of the DataPostprocessor-functionality of dealii, how would this 
> procedure know what the special calculation rule of the desired output 
> would look like?
>
> Concerning your second note, thanks for pointing me in the direction of  
> the Zienkiewicz-Zhu recovery procedure. The procedure is straightforward to 
> me, however its implementation in deal.ii not (yet). There are some 
> functions like the VectorTools::project() function, which is supposed to 
> give a L2-projections. 
> However, I am not quite sure how to use them. Therefore I tried to 
> implement the  Zienkiewicz-Zhu recovery procedure myself. But there is a 
> fundamental problem I could not solve so far. In elasticity problems for 2D 
> or 3D we have a vector-valued solution. Therefore, in 2D for a 
> triangulation with n Nodes we have 2n DoFs. If I want to add the projected 
> scalar function (e.g. the mises stress) in form of a command like
>
> data_out.add_data_vector(mises_stress_recovery, "mises_recovery");
>
> to the output, there is the problem that the mises_stress_recovery vector 
> only has a dimension of n, where as data_out wants a vector of dimension 
> 2n. How to deal with that issue?
> I tried to generate a new dof_handler via 
> DoFHandler dof_handler_recovery(this->triangulation);
> FESystem fe_recovery(FE_Q(this->fe_degree),1); // here dim=1 
> since we only have on dof per node
> MappingQ mapping_recovery(this->fe_degree);
> AffineConstraints constraints_recovery;
> dof_handler_recovery.distribute_dofs(fe_recovery);
>
> and then assemble a mass matrix and a rhs to solve for the 
> mises_stress_recovery vector. However, when looping over the cells, I am 
> not sure how to reference to the cells. Should I use the cells of the newly 
> created DoFHandler via 
> for (const auto  : dof_handler_recovery.active_cell_iterators())
> or via the DoFHandler of the original problem solution via
> for (const auto  : dof_handler.active_cell_iterators()).
>
> Do I really have to create a new DoFHandler? If not, how would one 
> assemble the mass matrix, since for that we need the an local_dof_indices 
> object like std::vector that points only to 
> every second dof index or so.
> I would be grateful for some tips or directions.
>
> 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/c0377cea-0a24-467f-bca9-c93be89dd8ben%40googlegroups.com.
#include 
#include 
#include 

#include 
#include 
#include 
#include 
#include 
#include 
#include 

#include 
#include 
#include 
#include 

#include 
#include 

#include 

#include 
#include 
#include 
#include 

#include 
#include 

#include 
#include 

double TOL = 10E-9;
namespace Step8{
  using namespace dealii;

  template 
  class ElasticProblem{
  public:
ElasticProblem();
void run();

  private:
double E=1.,nu=0.3;
void setup_system();
void assemble_system();
void solve();
void output_results() const;
Triangulation triangulation;
DoFHandlerdof_handler;
FESystem fe;
AffineConstraints constraints;
SparsityPattern  sparsity_pattern;
SparseMatrix system_matrix;
Vector solution;
Vector system_rhs;
  };

  template 
  void right_hand_

Re: [deal.II] Averaging derived solution quantities (mises stress)

[Konrad Schneider]

Dear Wolfgang,
yes the mesh is quite coarse, but I want to investigate the functionality 
of how deal.ii outputs stresses and how these stresses are then displayed 
via for instance paraview. To really see the effects I intentionally use 
the coarse mesh.
So far, I actually do create a (derived) DataPostprocessor-Object (see the 
thinned out and altered step-8 tutorial attached to this message) similar 
to the procedure in step-32. You wrote, that the gradients are computed at 
the vertices. How so? If I look at line 291 of my code or at step-32, there 
is a loop over all quadrature points (  for (unsigned int q = 0; q < 
n_quadrature_points; ++q) ) which gets called for every cell. In the scope 
of this loop the shape function gradients at the quadrature points are used 
to compute the desired output variable by a special calculation rule. 
However, these are not the vertices. If the gradients are somehow again 
computed at the vertices, which might be an internal procedure of the 
DataPostprocessor-functionality of dealii, how would this procedure know 
what the special calculation rule of the desired output would look like?

Concerning your second note, thanks for pointing me in the direction of  
the Zienkiewicz-Zhu recovery procedure. The procedure is straightforward to 
me, however its implementation in deal.ii not (yet). There are some 
functions like the VectorTools::project() function, which is supposed to 
give a L2-projections. 
However, I am not quite sure how to use them. Therefore I tried to 
implement the  Zienkiewicz-Zhu recovery procedure myself. But there is a 
fundamental problem I could not solve so far. In elasticity problems for 2D 
or 3D we have a vector-valued solution. Therefore, in 2D for a 
triangulation with n Nodes we have 2n DoFs. If I want to add the projected 
scalar function (e.g. the mises stress) in form of a command like

data_out.add_data_vector(mises_stress_recovery, "mises_recovery");

to the output, there is the problem that the mises_stress_recovery vector 
only has a dimension of n, where as data_out wants a vector of dimension 
2n. How to deal with that issue?
I tried to generate a new dof_handler via 
DoFHandler dof_handler_recovery(this->triangulation);
FESystem fe_recovery(FE_Q(this->fe_degree),1); // here dim=1 
since we only have on dof per node
MappingQ mapping_recovery(this->fe_degree);
AffineConstraints constraints_recovery;
dof_handler_recovery.distribute_dofs(fe_recovery);

and then assemble a mass matrix and a rhs to solve for the 
mises_stress_recovery vector. However, when looping over the cells, I am 
not sure how to reference to the cells. Should I use the cells of the newly 
created DoFHandler via 
for (const auto  : dof_handler_recovery.active_cell_iterators())
or via the DoFHandler of the original problem solution via
for (const auto  : dof_handler.active_cell_iterators()).

Do I really have to create a new DoFHandler? If not, how would one assemble 
the mass matrix, since for that we need the an local_dof_indices object 
like std::vector that points only to every second 
dof index or so.
I would be grateful for some tips or directions.

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/c711eb18-b84b-4120-a46e-5bfc76773ec0n%40googlegroups.com.
#include 
#include 
#include 
 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
 
#include 
#include 
#include 
#include 

#include 
#include 
 
#include 
 
#include 
#include 
#include 
#include 
 
#include 
#include 
 
#include 
#include 

double TOL = 10E-9;
namespace Step8{
using namespace dealii;
  
template 
class ElasticProblem{
  public:
ElasticProblem();
void run();
  
  private:
double E=1.,nu=0.3;
void setup_system();
void assemble_system();
void solve();
void output_results() const;
Triangulation triangulation;
DoFHandlerdof_handler;
FESystem fe;
AffineConstraints constraints;
SparsityPattern  sparsity_pattern;
SparseMatrix system_matrix;
Vector solution;
Vector system_rhs;
};
 
template 
void right_hand_side(const std::vector> ,std::vector> &  values){
  AssertDimension(values.size(), points.size());
  Assert(dim >= 2, ExcNotImplemented());
  Point point_1, point_2;
  point_1(0) = 0.5;
  point_2(0) = -0.5;
  for (unsigned int point_n = 0; point_n < points.size(); ++point_n)
{
  if (((points[point_n] - point_1).norm_square() < 0.2 * 0.2) ||
  ((points[point_n] - point_2).norm_square() < 0.2 * 0.2))
values[p

Re: [deal.II] output derived quanteties (gradients of vector valued solution like strain)

[Konrad Schneider]

Thanks. That was too easy. How should one know this...
Actually in the documentation 
(https://www.dealii.org/current/doxygen/deal.II/classDataPostprocessorTensor.html)
 
the command
 
data_out.write_vtk (output);

appears in the example. Maybe it should be changed to 

data_out.write_vtu (output);

Furthermore, may I suggest, that it might be of some value to add the data 
output of stresses to one of the elasticity tutorials (maybe add it to 
step-8), since evaluating stresses (especially von mises stresses) is one 
of the major tasks in "pure" solid mechanics plus there is this pitfall 
that we have to use the vtu-format for tensor data.

But thanks anyway for the quick response (again).

Best regards
Konrad

Wolfgang Bangerth schrieb am Freitag, 10. März 2023 um 20:35:09 UTC+1:

> On 3/10/23 07:41, Konrad Schneider wrote:
> > 
> > I am using deal.ii version 9.4 and am wondering what I did wrong here. 
> Do you 
> > have any suggestions?
>
> The VTK writer simply doesn't support writing tensor data. Use the VTU 
> writer 
> instead.
>
> Best
> W.
>
> -- 
> 
> Wolfgang Bangerth email: bang...@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/c5d4f3f2-e741-4008-a712-e7d108fb7064n%40googlegroups.com.

[deal.II] output derived quanteties (gradients of vector valued solution like strain)

[Konrad Schneider]

Dear Forum,
I tried to follow the instructions of the DataPostprocessorTensor class to 
output the gradient of a vector valued solution. Therefore I tried to adopt 
step-8. However, I am getting the following error message:

An error occurred in line <5640> of file 
 in function
void dealii::DataOutBase::write_vtk(const 
std::vector >&, const 
std::vector >&, const 
std::vector, 
std::allocator >, 
dealii::DataComponentInterpretation::DataComponentInterpretation> >&, const 
dealii::DataOutBase::VtkFlags&, std::ostream&) [with int dim = 2; int 
spacedim = 2; std::ostream = std::basic_ostream]
The violated condition was: 
std::get<3>(nonscalar_data_range) != 
DataComponentInterpretation::component_is_part_of_tensor
Additional information: 
You are trying to use functionality in deal.II that is currently not
implemented. In many cases, this indicates that there simply didn't
appear much of a need for it, or that the author of the original code
did not have the time to implement a particular case. If you hit this
exception, it is therefore worth the time to look into the code to
find out whether you may be able to implement the missing
functionality. If you do, please consider providing a patch to the
deal.II development sources (see the deal.II website on how to
contribute).

I am using deal.ii version 9.4 and am wondering what I did wrong here. Do 
you have any suggestions?
Best regards 
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/a04ecfe4-3518-4253-930e-addb6f619fd4n%40googlegroups.com.
#include 
#include 
#include 
 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
 
#include 
#include 
#include 
 
#include 
#include 
 
#include 
 
#include 
#include 
#include 
#include 
 
#include 
#include 
 
#include 
#include 
 
namespace Step8
{
  using namespace dealii;
  
template 
class ElasticProblem{
  public:
ElasticProblem();
void run();
  
  private:
void setup_system();
void assemble_system();
void solve();
void refine_grid();
void output_results() const;
Triangulation triangulation;
DoFHandlerdof_handler;
FESystem fe;
AffineConstraints constraints;
SparsityPattern  sparsity_pattern;
SparseMatrix system_matrix;
Vector solution;
Vector system_rhs;
};
 
template 
void right_hand_side(const std::vector> ,std::vector> &  values){
  AssertDimension(values.size(), points.size());
  Assert(dim >= 2, ExcNotImplemented());
  Point point_1, point_2;
  point_1(0) = 0.5;
  point_2(0) = -0.5;
  for (unsigned int point_n = 0; point_n < points.size(); ++point_n)
{
  if (((points[point_n] - point_1).norm_square() < 0.2 * 0.2) ||
  ((points[point_n] - point_2).norm_square() < 0.2 * 0.2))
values[point_n][0] = 1.0;
  else
values[point_n][0] = 0.0;

  if (points[point_n].norm_square() < 0.2 * 0.2)
values[point_n][1] = 1.0;
  else
values[point_n][1] = 0.0;
}
}

template 
ElasticProblem::ElasticProblem()
  : dof_handler(triangulation)
  , fe(FE_Q(1), dim)
{}
 
template 
void ElasticProblem::setup_system()
{
  dof_handler.distribute_dofs(fe);
  solution.reinit(dof_handler.n_dofs());
  system_rhs.reinit(dof_handler.n_dofs());

  constraints.clear();
  DoFTools::make_hanging_node_constraints(dof_handler, constraints);
  VectorTools::interpolate_boundary_values(dof_handler,
0,
Functions::ZeroFunction(dim),
constraints);
  constraints.close();

  DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
  DoFTools::make_sparsity_pattern(dof_handler,
  dsp,
  constraints,
  /*keep_constrained_dofs = */ false);
  sparsity_pattern.copy_from(dsp);

  system_matrix.reinit(sparsity_pattern);
}
 
template 
void ElasticProblem::assemble_system()
{
  QGauss quadrature_formula(fe.degree + 1);
  FEValues fe_values(fe,
  quadrature_formula,
  update_values | update_gradients |
  update_quadrature_points | update_JxW_values);
  const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
  const unsigned int n_q_points= quadrature_formula.size();

  FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell);
  Vector cell_rhs(dofs_per_cell);

  std

Re: [deal.II] Required mapping member initialization of type MappingQ

[Konrad Schneider]

Dear Wolfgang,

thanks for clearing this up and responding so fast. Your explanation makes 
sense to me.

Best
Konrad

Wolfgang Bangerth schrieb am Mittwoch, 8. März 2023 um 15:03:13 UTC+1:

> On 3/8/23 06:34, Konrad Schneider wrote:
> > 
> > Only if I uncomment line 30 of my code (// , mapping(mapping_degree)) 
> the 
> > program compiles. Why is that? In my understanding, the program should 
> compile 
> > just fine without the initialization of the mapping member variable 
> since 
> > there is no requirement of initializing member variables upon 
> constructing an 
> > object, is would be the case for my member variable unsigned int test.
>
> This isn't true. In a constructor, *all* member variables [1] are 
> initialized 
> by calling their constructors. If you explicitly list a member in the 
> initializer list (after the ':'), then the constructor arguments so 
> specified 
> are taken. If you don't explicitly list a variable, then the default 
> constructor of the variable's class is called. But the latter only works 
> if 
> the class *has* a default constructor. It turns out that MappingQ does 
> not: 
> The only constructor there is requires you to provide the polynomial 
> degree. 
> As a consequence, the compiler cannot do a default-initialization, and you 
> are 
> forced to explicitly list the variable in the initializer list.
>
> Best
> W.
>
>
> [1] This is not quite true: all *class type* member variables are 
> initialized, 
> whereas variables of built-in types like int, double, etc are not unless 
> you 
> explicitly initialize them.
>
> -- 
> 
> Wolfgang Bangerth email: bang...@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/8a6d8887-dd98-4945-9325-f51731f10e3dn%40googlegroups.com.

[deal.II] Required mapping member initialization of type MappingQ

[Konrad Schneider]

Dear all,
I am fairly new to deal.ii and have a question concerning the 
MappingQ-class of deal.ii
Why do I have to initialize the private mapping member of type 
MappingQ to get a compilable code? To illustrated what I mean, I did 
thin out the step-11 tutorial in the following way:

 
#include 
#include 
#include 
#include 
#include 
#include 
using namespace dealii;

template 
class LaplaceProblem{
public:
LaplaceProblem(const unsigned int mapping_degree,
const unsigned int fe_degree=1);
void run();
private:
Triangulation triangulation;
const unsigned int fe_degree;
FE_Q fe;
DoFHandler dof_handler;
MappingQ mapping;
unsigned int test;
};

template 
LaplaceProblem::LaplaceProblem(const unsigned int mapping_degree,
const unsigned int fe_degree)
: fe_degree(fe_degree)
, fe(fe_degree)
, dof_handler(triangulation)
// , mapping(mapping_degree)
{
}

template 
void LaplaceProblem::run(){
}
int main(){
const unsigned int mapping_degree=2;
LaplaceProblem<2>(mapping_degree).run();
}


Only if I uncomment line 30 of my code (// , mapping(mapping_degree)) the 
program compiles. Why is that? In my understanding, the program should 
compile just fine without the initialization of the mapping member variable 
since there is no requirement of initializing member variables upon 
constructing an object, is would be the case for my member variable unsigned 
int test  .
I would be grateful for any explanation.
Many thanks in advance.

-- 
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/92911e1a-8ca1-425f-a859-9b5a6e877b4dn%40googlegroups.com.

[deal.II] MPI: Evaluating several distributed functions on each node

[Konrad Simon]

Dear all,

I am facing a little challenge: Suppose I have several objects myFunction 
of a certain class MyFunction derived from a Function object on 
each process in my MPI universe. The number of MyFunction objects may vary 
from process to process.
I am trying to evaluate each object at certain points p that live on a 
distributed::Triangualtion. The evaluation points are different for 
each MyFunction object. A point may be contained in a cell owned by another 
process whose id is known, so I know where to send it to to ask a 
MyFunction object on that specific process to evaluate it for me and send 
the value back.

Now my question: How do I deal with several MyFunction objects on one 
process (I can identify these by a hash). Is it a good idea to have a class 
with a set of static variables and methods on each node to manage the 
communication and distribute necessary information within the process using 
the specific MyFunction hash? Does that make sense?

Let me stress again, the difficulty here is that I have several objects 
that may ask for values one by one or even in a (shared memory) parallel 
manner that may or may not trigger communication (sort of some_to_some 
style as implemented in Deal.II). 

Thanks in advance :-)

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/4fb6b294-aa7f-4dc2-b823-e42247bf2302n%40googlegroups.com.

[deal.II] Triangulation from lower-dimensional object

[Konrad Simon]

Dear all,

I was recently looking for a solution to the following problem:

I have a triangulation, e.g., of a unit cube. I need to solve a 
lower-dimensional (FEM-like) problem on all faces with BCs on edges and for 
each face I need to solve a problem on all edges with BCs at vertices.

Can I somehow extract a Triangulation form a Triangulation by 
providing, for example, the boundary id on the Triangulation level?

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/ae0eae00-a572-4036-b20e-31f45dd70892n%40googlegroups.com.

[deal.II] Elasticity Preconditioner

[Konrad Simon]

Dear all, 

A student and me are trying to deal with (parallel, static) linearized 
elasticity similar to step-8.
However, the problem is slightly different since the Lame parameters are 
functions instead of constants. 
Can anyone recommend good solvers & preconditioners for the resulting 
system? Our boundary conditions fix all degrees of freedom in the kernel.

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/eba0cc0f-6028-4fde-9cee-307e09f8d2ban%40googlegroups.com.

[deal.II] Re: Question about deal.II flexibility

[Konrad Simon]

Dear Abdo,

I do think that Deal.II can help you with solving your problem. If you are 
just starting with Deal.II have a look at the tutorials 
<https://dealii.org/developer/doxygen/deal.II/Tutorial.html>. It is always 
a good idea to start with some code basis that other people have shared and 
that was thoroughly tested. As to your future problem, please have a look 
at the code gallery 
<https://dealii.org/developer/doxygen/deal.II/CodeGallery.html> as well. 
There are codes that other users shared dealing with visco-plastic 
materials, in particular (but not only Jean-Paul).

Cheers,
Konrad

On Saturday, June 19, 2021 at 10:23:37 PM UTC+2 aael...@ncsu.edu wrote:

> Hi
> I am Abdo, Ph.D. student at NC state university and was looking for a FE 
> package library that is flexible and can be used in modeling and analyzing 
> 2D and 3D elastic wave equation propagation within an incompressible 
> medium. Also, at some point, I might need to add the viscoelastic behavior 
> to the model. So, I am not sure if deal.II is flexible enough to add, 
> modify some feature in it or not such that I can do what I have mentioned 
> earlier. So, Can you please give me some pointers or advice regarding 
> this?, your help is much appreciated.
>
> Best regards.
> Abdo.
>

-- 
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/481afb9d-fa92-4eae-b369-0881cbd7d70fn%40googlegroups.com.

Re: [deal.II] p4est user pointer and distributed::triangulation

[Konrad Simon]

Thank you, Wolfgang! This is exactly what I am looking for. I don‘t need to 
access the pointer from outside, I am trying to extend the p4est-interface 
itself.

Best,
Konrad

On Friday, May 14, 2021 at 12:42:48 AM UTC+2 Wolfgang Bangerth wrote:

> On 5/13/21 3:05 PM, Konrad Simon wrote:
> > 
> > I am currently writing an interface to some p4est functions. The 
> structure of 
> > p4est makes it often necessary to pass data around through a user 
> pointer of 
> > whatever type (void*). When creating a new triangulation this pointer is 
> set 
> > to „this“ (the triangulation itself), see for example
> > 
> > 
> https://www.dealii.org/current/doxygen/deal.II/distributed_2tria_8cc_source.html#l2989
>  
> > <
> https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.dealii.org%2Fcurrent%2Fdoxygen%2Fdeal.II%2Fdistributed_2tria_8cc_source.html%23l2989=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cb4cadd28da84445fd29308d91652e2c3%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637565367510208994%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000=v2G100jGfI6SbGwV35raSFjMxfnJMJurP2F2qAc61Ls%3D=0
> >
> > 
> > Is this pointer used in other interfaces somehow? p4est itself does not 
> touch 
> > it so I am wondering if I can reset it to anything I‘d like.
>
> Konrad,
> we use this user pointer in many of the functions that are called back 
> from 
> p4est. Take a look at
> RefineAndCoarsenList::refine_callback
> for example (line 761 of the file), and many other functions that you can 
> find 
> by searching for
> forest->user_pointer
> in this file.
>
> We consider the p4est forest object as internal to the p::d::T class, so I 
> would expect that you can't access it, or expect that you can use any of 
> its 
> content for anything other than what the p::d::T class does with it.
>
> Best
> W.
>
>
> -- 
> 
> Wolfgang Bangerth email: bang...@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/7fbd1adf-c6d4-4656-ab6e-7c833e5e92cen%40googlegroups.com.

[deal.II] p4est user pointer and distributed::triangulation

[Konrad Simon]

Dear all, 

I am currently writing an interface to some p4est functions. The structure 
of p4est makes it often necessary to pass data around through a user 
pointer of whatever type (void*). When creating a new triangulation this 
pointer is set to „this“ (the triangulation itself), see for example

https://www.dealii.org/current/doxygen/deal.II/distributed_2tria_8cc_source.html#l2989

Is this pointer used in other interfaces somehow? p4est itself does not 
touch it so I am wondering if I can reset it to anything I‘d like.

Any Ideas?

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/f23950cf-7688-45d2-b986-140ccc3a58f9n%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 
<https://dealii.org/current/doxygen/deal.II/group__constraints.html#gacdba6b0a8a06f9c9ceff9f11ef05a203>
 
or project_boundary_values_div_conforming 
<https://dealii.org/current/doxygen/deal.II/group__constraints.html#gafb7e51448054f5a4cdc95c7e39dfc0e0>

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 
<https://epubs.siam.org/doi/book/10.1137/1.9781611975543>. 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.

Re: [deal.II] interpolate FE_Nelelec

[Konrad Simon]

Dear John,

As Jean-Paul pointed out the entries of the solution vector for a Nedelec 
element have a different meaning. The entries are weights that result from 
edge moments (and in case of higher order also face moments and volume 
moments), i.e, integrals on edges. Nedelec elements have a deep geometric 
meaning: they discretize quantities that you can integrate (a physicist 
would say "measure") along 1D-lines in 3D.

Btw, if you are using Nedelec in 2d be aware that it has some difficulties 
on complicated meshes. In 3D the lowest order is fine. NedelecSZ (for 2D 
lowest order and 3D all orders), as Jean-Paul pointed out, are another 
option. They are meant to by-pass certain mesh orientation issues on 
complicated geometries (I can tell you a bit about that since currently I 
am chasing some problems there).

Best,
Konrad

On Tuesday, April 20, 2021 at 9:23:23 PM UTC+2 Jean-Paul Pelteret wrote:

> Hi John,
>
> Unfortunately, you’ve fallen into the trap of confusing what the entries 
> in the solution vector mean for the different element types. For Nedelec 
> elements, they really are coefficients of a polynomial function, so you 
> can’t simply set each coefficient to 1 to visualise the shape function 
> (like one might be inclined to try for Lagrange type polynomials). If you 
> want a little piece of code that will plot the Nedelec shape functions for 
> you, then take a look over here:
> https://github.com/dealii/dealii/issues/8634#issuecomment-595715677
>
> BTW. In case you’re not aware, a newer (and probably more robust) 
> implementation of a Nedelec element is the FE_NedelecSZ 
> <https://dealii.org/current/doxygen/deal.II/classFE__NedelecSZ.html> that 
> uses hierarchical shape functions to form the higher order bases. The PhD 
> thesis of S. Zaglmayr that describes the element is mentioned in the class 
> documentation.  
>
> Best,
> Jean-Paul
>
> On 20. Apr 2021, at 15:47, John Smith  wrote:
>
>
> Dear Konrad,
>
> Thank you for your help! I have recently modified the step-3 tutorial 
> program to save FE_Nedelec<2> fe(0) shape function into *.csv files. It is 
> attached to this message. The result of this program is not what I would 
> expect it to be. I have also attached a pdf file which compares my 
> expectations (marked as “Expectations” in the pdf file) with the results of 
> the modified step-3 program (marked as “Reality” in the pdf file). My 
> expectations are based on the literature I have. My question is: why the 
> shape functions of FE_Nedelec<2> fe(0) differ from my expectations? I guess 
> my expectations are all wrong. I would also appreciate if you will share 
> with me a reference to a paper or a book, so I can extend my modest library.
>
> Thanks!
>
> John
>
> P.S. I have sent you an e-mail as you have requested. Could you please 
> check your spam folder? 
>
>
>
> On Monday, April 19, 2021 at 8:47:29 AM UTC+2 Konrad Simon wrote:
>
>> Dear John,
>>
>> I had a look at the pdf you sent. I noticed some conceptual 
>> inconsistencies that are important for the discretization. Let me start 
>> like this: The curl-curl-problem that you are trying to solve is - 
>> according to your description of the gauge - actually a curl-curl + 
>> grad-div problem and hence a Laplace-problem that describes a 
>> Helmholtz-decomposition. In other words in order to get a well-posed 
>> simulation problem you need more boundary conditions. You already noticed 
>> that since your curl-curl-matrix is singular (the kernel is quite 
>> complicated and contains all longitudinal waves). 
>>
>> You have of course the option to solve (0.0.2) with a Krylov-solver but 
>> then you need to make sure that the right-hand side is in the orthogonal 
>> complement of this kernel before each iteration which is quite difficult. I 
>> would not recommend that option.
>>
>> Another point is that if you have a solution for your field A like in 
>> (0.0.4) you can not have a similar representation for T. This is 
>> mathematically not possible.
>>
>> Since you not care about the gauge let me tell you how I would tackle 
>> this: The reason  is  - to come back to my original point - you are missing 
>> a boundary condition that makes you gauge unique. Since you apply natural 
>> boundary conditions (BCs) on A you must do so as well to determine div(A). 
>> This second BC for A is applied to a different part of the system that is 
>> usually neglected in the literature and A is partly determined from this 
>> part (this part then describes the missing transversal waves). The 
>> conditions you mention on curl(A) are some what contradictory to the space 
>> in whic

Re: [deal.II] Re: interpolate FE_Nelelec

[Konrad Simon]

Dear John,

I had a look at the pdf you sent. I noticed some conceptual inconsistencies 
that are important for the discretization. Let me start like this: The 
curl-curl-problem that you are trying to solve is - according to your 
description of the gauge - actually a curl-curl + grad-div problem and 
hence a Laplace-problem that describes a Helmholtz-decomposition. In other 
words in order to get a well-posed simulation problem you need more 
boundary conditions. You already noticed that since your curl-curl-matrix 
is singular (the kernel is quite complicated and contains all longitudinal 
waves). 

You have of course the option to solve (0.0.2) with a Krylov-solver but 
then you need to make sure that the right-hand side is in the orthogonal 
complement of this kernel before each iteration which is quite difficult. I 
would not recommend that option.

Another point is that if you have a solution for your field A like in 
(0.0.4) you can not have a similar representation for T. This is 
mathematically not possible.

Since you not care about the gauge let me tell you how I would tackle this: 
The reason  is  - to come back to my original point - you are missing a 
boundary condition that makes you gauge unique. Since you apply natural 
boundary conditions (BCs) on A you must do so as well to determine div(A). 
This second BC for A is applied to a different part of the system that is 
usually neglected in the literature and A is partly determined from this 
part (this part then describes the missing transversal waves). The 
conditions you mention on curl(A) are some what contradictory to the space 
in which A lives. Nedelec elements which you must use (you can not use FE_Q 
to enforce the conditions) cannot generate your desired T_0.

There are some principles when discretizing these problems (which are not 
obvious) that you MUST adhere to (choice of finite elements, boundary 
conditions, what is the exact system etc) if you want to get a stable 
solution and these are only understood recently. I am solving very similar 
problems (with Deal.II) in a fluid context and will be happy to share my 
experiences with you. Just email me: konrad.si...@uni-hamburg.de

Best,
Konrad

On Monday, April 19, 2021 at 5:51:17 AM UTC+2 Wolfgang Bangerth wrote:

> On 4/18/21 12:24 PM, John Smith wrote:
> > 
> > Thank you for your reply. It is a bit difficult to read formulas in 
> text. So I 
> > have put a few questions I have in a pdf file. Formulas are better 
> there. It 
> > is attached to this message. May I ask you to have a look at it?
>
> John:
>
> If I understand your first question right, then you are given a vector 
> field J 
> and you are looking for a vector field T so that
> curl T = J
> I don't really have anything to offer to this kind of question. There are 
> many 
> vector fields T that satisfy this because of the large null space of curl. 
> You 
> have a number of condition you would like to "approximate" but there are 
> many 
> ways to "approximate" something. In essence, you have two goals: To 
> satisfy 
> the equation above and to approximate something. You have to come up with 
> a 
> way to weigh these two goals. For example, you could look to minimize
> F(T) = ||curl T - J||^2 + alpha ||T-T_desired||^2
> where T_desired is the right hand side of (0.0.5) and you have to 
> determine 
> what alpha should be.
>
> As for your other questions: (0.0.9) is indeed called "interpolation"
>
> The issue with the Nedelec element (as opposed to a Q(k)^dim) field is 
> that 
> for the Nedelec element,
> \vec phi(x_j)
> is not independent of
> \vec phi(x_k)
> and so you can't choose the matrix in (0.0.8) as the identity matrix. You 
> realize this in (0.0.13). The point I was trying to make is that (0.0.13) 
> cannot be exactly satisfied if T(x_j) is not a vector parallel to \hat 
> e_j, 
> which in general it will not be. You have to come up with a different 
> notion 
> of what it means to solve (0.0.13) because you cannot expect the left and 
> right hand side to be equal.
>
> Best
> W>
>
> -- 
> 
> Wolfgang Bangerth email: bang...@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/75104172-4893-4ae6-af8e-8ec2c12d335an%40googlegroups.com.

Re: [deal.II] 2d meshes and orientation

[Konrad Simon]

Dear Wolfgang,

On Tuesday, April 6, 2021 at 5:53:15 AM UTC+2 Wolfgang Bangerth wrote:

> It's possible this is indeed a bug. In most cases, we run the 2d meshes 
> through the mesh orienter, so we rarely see 2d meshes that are not 
> correctly 
> oriented and it wouldn't surprise me if we have the kind of bug you 
> describe 
> where things almost but not quite work correctly.


I see. Many thanks for the input. This is my impression as well. These 
things rarely happen in 2d.
One example would be simple periodic meshes which cause line orientations 
not to match (this case is caught). In more
complicated scenarios involving periodicity and rotations of faces it can 
happen though that lines do match but normals do not.


> The questions then are (i) where does this need to be fixed, and (ii) how 
> many 
> tests fail if we apply the fix? One of the things I learned (and I think 
> you 
> did as well) is that for bugs in these kinds of layers, it is quite 
> possible 
> that one ends up with 100 failing tests and then it takes a *lot* of work 
> to 
> look through all of these to figure out whether the new output is correct 
> or 
> not. But then, as mentioned above, it's also possible that only one or two 
> tests fail because we almost never see these kinds of meshes.
>

Sound worth trying. Please see also this new PR 
 which is work in progress.

-- 
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/0125f6dc-b61c-4db6-84da-b101b6759f1bn%40googlegroups.com.

[deal.II] 2d meshes and orientation

[Konrad Simon]

Dear deal.ii community,

I stumbled over something in 2d meshes that confuses me a bit. Supoose I 
have 2 squares (a left and a right one) sharing one edge. If everything is 
fine all lines and faces have standard orientation, i.e., 
line_orientation==true and face_orientation==true for all lines and faces. 
Note that lines and faces are NOT the same geometric entity here since 
their orientation is determined differently (faces have normals and lines 
have tangents).

Now, if I rotate the right square by 90 degrees I get (on the shared edge) 
non-matching face normals but all line orientations do match.

If I rotate the right square by 180 degrees I get (on the shared edge) 
non-matching face normals and non-matching line orientations.

If I rotate the right square by 270 degrees I get (on the shared edge) 
matching face normals and non-matching line orientations.

This is at least my intuition. When I query these information for each 
rotation case I see that my intuition is confirmed for lines but the face 
orientations are always true.

Is this a bug?

I am asking because it seems that 2d RaviartThomas elements (which make use 
of face normals) are broken for some edge cases (3d RT I could fix).

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/f739c450-c37d-4d24-bf5b-2778dec5d501n%40googlegroups.com.

[deal.II] Re: Need help installing deal.ii with p4est and Trilinos

[Konrad Simon]

Dear Budhyant,, 

Hard, to say what should be the paths on your system. Dealing with all 
these large dependencies can be quite cumbersome. I would say the way you 
install Deal.II depends on what you want to do with it. Let me tell you how 
I would approach this.

1. If you just want to develop a parallel application (using all the 
packages that you mention) on your local laptop it is in my opinion usually 
enough to use the packages that are provided by the OS package manager. 
Simply use apt install if it is your own machine and you have super user 
permissions. Then cmake should find these packages in the PATH and you just 
need to tell cmake to use them through, e.g., DEAL_II_WITH_P4EST = ON etc, 
for the build of Deal.II. This is usually sufficient for development since 
for running large 3D problems a single machine is anyway not enough. So no 
need to optimize a build for speed

2. If you do not have root permissions (or if you want to do some larger 
computations on a cluster, or you don't want to change the system packages) 
installing Deal.II using spack <https://spack.io/> is quite a nice option. 
Spack will resolve all dependencies, optimize the builds for your system 
and can be customized nicely (you can also link against MPI libraries and 
other stuff provided by the target system's admin). This is the option that 
I use.

Hope that helps.
Best regards,
Konrad

On Friday, February 19, 2021 at 10:07:00 AM UTC+1 venepalli...@gmail.com 
wrote:

> Dear All:
>
> I would appreciate it if you could please help me with the installation on 
> Ubuntu 20. 
>
> I have installed p4est, petsc, Trilinos and metis in "~/src/", they are 
> working individually. 
>
> Cmake succeeds in recognizing petsc and metis, however, it fails to 
> recognize the p4est and Trilinos installations. 
>
> *cmake -DDEAL_II_WITH_MPI=ON -DP4EST_DIR=~/src/p4est/build-mpi/ 
> -DPETSC_DIR=~/src/petsc-3.14.4/ -DPETSC_ARCH=amd64 
> -DTRILINOS_DIR=~/src/Trilinos/MY_BUILD/ -DMETIS_DIR=~/src/metis-5.1.0/ 
> -DCMAKE_INSTALL_PREFIX=~/src/dealii-9.2.0/ ..*
>
> What should the path be in my case?
>
> Thank you.
>
> Regards,
> Budhyant Venepalli
>

-- 
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/af5964fe-cfb1-4be7-a679-0235cdbf2b08n%40googlegroups.com.

Re: [deal.II] Face and line orientation dependent dof permutations in FE_Q

[Konrad Simon]

> 0 1 3 4 (thinks the edge is 1->4) 
> 5 4 2 1 (thinks the edge is 4->1) 
> or any other rotation. A good test would probably just try all 2x4 
> possible 
> rotations of both cells. 
>

Yes, I submitted another pull request containing exactly this.
 

>
> 2/ How do you check whether things are right or not? There are many 
> aspects, 
> but one I've recently thought would be useful to try out is 
> FEInterfaceValues. 
> Try this: assign a DoFHandler with a continuous but possibly high-order 
> element with multiple DoFs per edge. Create a solution vector and give its 
> elements random elements. This should correspond to some finite element 
> field, 
> and because it is associated with a continuous FE, the jump across the 
> common 
> edge should be zero. 
>
> So create an FEInterfaceValues object and check that indeed the jump of 
> the 
> solution (or, more simply, of each shape function) at every quadrature 
> point 
> along that face is really zero. If it is not, there's a bug somewhere. 
>

Thanks for the hint with  FEInterfaceValues. This I did not think about but 
probably this is a good think for the test suite. 
I used visual debugging, a modified step-20 on problematic meshes and 
applications with a known solution. I could reproduce
expected convergence rates for RaviartThomas.
 

>
> My best guess is that with such small test programs, you'd uncover cases 
> that 
> don't work and that could, because they don't involve solving whole 
> problems 
> or more complicated set ups, be easier to debug and/or fix. Everything 
> that 
> turns out to work as expected should then be put into the test suite as 
> soon 
> as you've verified that it works -- we accept small additions to the test 
> suite pretty regularly, so feel free to open lots of pull requests in this 
> regard. 
>

This is how I fixed RaviartThomas. Works nicely. 

>
> As mentioned above, regrettably, we've lost collective knowledge about 
> these 
> issues over the years, and that translates into long response times :-( 
> I'm 
> sorry you get caught up in this. I'd love to have fixed this long ago, and 
> would have loved to help you more with it as well, but time availability 
> is 
> not always on my side :-(( 
>

I actually need all these vector elements so we both have an interest to 
fix it. 
Looking at the current implementation it seems like the developer who did 
it had good ideas but left the
work in the middle (there are many TODO comments and missing 
implementations).
There is a lot of work to do. Maybe I can help a bit here.

Thank you a lot, Wolfgang, for your elaborate answer :-)

Cheers,
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/652cf426-6b5d-4f3f-bfce-f1266e171d20n%40googlegroups.com.

[deal.II] Face and line orientation dependent dof permutations in FE_Q

[Konrad Simon]

Dear Deal.II community,

I am currently fixing an orientation issue for some (non-primitive) vector 
elements with the help of two other Deal.II developers/maintainers.

I have a quick question: On complicated meshes some faces can have a 
non-standard orientation or they can be rotated against the neighboring 
face. This causes inconsistencies for some vector elements.

In principle this potentially causes inconsistencies for FE_Q elements as 
well (for all dofs located at vertices, lines and faces). But that is being 
taken care of in the library, so everything seems fine.

In order to understand the issue better for other elements: where in the 
library is this being taken care of? I have difficulties to find that. 

Remark: I found a table in (fe.h/fe.cc) that derived elements need to 
implement that takes care of local dof permutations on quads (not full 
faces and hence not line dofs or vertex dofs). Now, lines can have 
different orientations but lines also permute within a cell upon face 
orientation/rotation changes.

Can anyone point me to the place in Deal.II where this is being taken care 
of for FE_Q?

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/7c3e0787-1b1c-465e-9747-3f82ff677d95n%40googlegroups.com.

Re: [deal.II] Mean value of component of TrilinosWrappers::MPI::BlockVector

[Konrad Simon]

>
> Then you need to create a vector with locally relevant entries as ghosts 
> and 
> copy your fully distributed vector to it. There is no other way if you 
> want to 
> use that function. 
>
> But if your goal is to fix the pressure in a Stokes problem, it doesn't 
> have 
> to be the integral mean of the pressure that is zero. It could just be the 
> arithmetic mean, which you can compute without access to ghost elements.  
>

Many  thanks, Wolfgang, that works! :-)

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/9849962a-5b13-4fb3-bcb5-bc74e045f2a5n%40googlegroups.com.

Re: [deal.II] Re: Mean value of component of TrilinosWrappers::MPI::BlockVector

[Konrad Simon]

> That looks like an unrelated error. Can you create a small testcase for 
> this 
> issue here? 
>

I will try to come up with an example. 

Thank you again and best regards,
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/83b966c1-4e3a-477c-827a-9d17f13682d8n%40googlegroups.com.

Re: [deal.II] Mean value of component of TrilinosWrappers::MPI::BlockVector

[Konrad Simon]

>
> I suspect you are passing a fully distributed vector to that function, but 
> it 
> needs to read ghost elements of the vector. Have you tried copying the 
> vector 
> into a locally_relevant vector, and passing that to the function in 
> question? 
>

Thank you, Wolfgang, that was the issue.  I am using the mean value 
function at several different places in my code. One of them is inside a 
preconditioner vmult that accepts only a fully distributed vector. Is there 
also a simple workaround? I have a patch since I only need the pressure 
which is the last component of the distributed vector in each locally owned 
cell but the patch is ugly...

-- 
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/1d0a185d-84af-4be5-8d2c-9b7979528c1cn%40googlegroups.com.

[deal.II] Re: Mean value of component of TrilinosWrappers::MPI::BlockVector

[Konrad Simon]

I am pretty sure that this is related to the following problem:

When I try to run AMR with this FESystem (Nedelec-RaviartThomas-DGQ) with 
more than one MPI rank I get this error

"An error occurred in line <1626> of file 
 
in function 
   void 
dealii::AffineConstraints::add_line(dealii::AffineConstraints::size_type)
 
[with number = double; dealii::AffineConstraints::size_type
= unsigned int] 
The violated condition was:  
   line_n != numbers::invalid_size_type 
Additional information:  
   This exception -- which is used in many places in the library -- usually 
indicates that some condition which the author of the code thought must be 
satisfied at a
certain point in an algorithm, is not fulfilled. An example would be that 
the first part of an algorithm sorts elements of an array in ascending 
order, and a second 
part of the algorithm later encounters an element that is not larger than 
the previous one. 

There is usually not very much you can do if you encounter such an 
exception since it indicates an error in deal.II, not in your own program. 
Try to come up with the 
smallest possible program that still demonstrates the error and contact the 
deal.II mailing lists with it to obtain help.
"

which is thrown in AffineConstraints::add_line 
<https://dealii.org/current/doxygen/deal.II/affine__constraints_8h_source.html#l01618>
 
(see the second assert in line 1826) after being called by 
make_hanging_node_constraints. It seems that deal.ii is trying to setup 
constraints on dofs that are not owned by the MPI process.

I read a little bit in the discussions referring to Nedelec and refinement 
- seems a difficult issue. Did anyone run into this, too?

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/a3b92b45-1767-42ee-87ce-3e668a331a93n%40googlegroups.com.

[deal.II] Mean value of component of TrilinosWrappers::MPI::BlockVector

[Konrad Simon]

Dear all,

I am trying to compute the mean value of the pressure component of a 
Trilinos block vector with three blocks and 7 components 
(vorticity-velocity-pressure).

Using one MPI rank is fine but if I use more I get the error

"An error occurred in line <666> of file 
 in 
function 
   dealii::TrilinosScalar 
dealii::TrilinosWrappers::MPI::Vector::operator()(dealii::TrilinosWrappers::MPI::Vector::size_type)
 
const 
The violated condition was:  
   false 
Additional information:  
   You tried to access element 3023 of a distributed vector, but this 
element is not stored on the current processor. Note: There are 4456 
elements stored on the current processor from within the range 5544 through 
 but Trilinos vectors need not store contiguous ranges on each 
processor, and not every element in this range may in fact be stored 
locally."

Now, I know what that means. Is there anything I can do? A workaround?

Thanks in advance and best regards,
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/15e015b4-af0f-466e-b628-f3777ab64327n%40googlegroups.com.

Re: [deal.II] Div-conforming elements in 3D

[Konrad Simon]

Dear Jean-Paul,

Many thanks for your reply.

On Tuesday, December 29, 2020 at 9:31:49 PM UTC+1 Jean-Paul Pelteret wrote:

> Hi Konrad,
>
> I'm sorry for taking some time to reply. To be honest, the inner working 
> of the FE classes is not something that I've ever had the time or 
> opportunity to dig into, so I'm really not well positioned to answer many 
> of your questions. I'm glad that you've submitted a PR for the work that 
> you've done thus far, and I hope that you can be guided by someone who has 
> more familiarity with this part of the library.
>
 
It will be my pleasure if other people will be able to benefit from 
implementations that I had to add and to contribute it as well as I 
benefited from other peoples‘ work. Let‘s hope that the stuff I asked to 
add makes sense. 
 

>
> What I can say, though, is that the various FE classes have been 
> implemented by numerous people over time, and the logic that dictates the 
> way that the DoFs are assigned varies between each FE. We have always 
> considered this to be perfectly fine -- this is, to the average person, an 
> implementational detail, and the interface that we have to interrogate DoFs 
> (e.g. which vector component they're associated with, whether or not they 
> have support points, etc...) seems to be sufficiently rich to most users. 
> So, as the end user, you shouldn't be writing code that would care that the 
> FE_Q element orders DoFs like this: [N_{0,x} N_{0,y} N_{0,z} N_{1,x} ...] ; 
> while IIRC the FE_DGQ element orders DoFs like this: [N_{0,x} N_{1,x} 
> N_{2,x} ... N_{0,y} N_{1,y} ...] .
>
> So, returning to your original set of questions: In general, I don't 
> believe that it can be assumed that all cell vertex DoFs are enumerated 
> before line/edges DoFs, before face DoFs, before cell interior DoFs. 
> Through the FiniteElementData class, the FiniteElement base class seems to 
> have a set of n_dofs_per_[vertex,line,face,cell]_dofs() 
> <https://dealii.org/developer/doxygen/deal.II/classFiniteElementData.html#a0d1f06778f87a7606cc2e0d5338b41ab>
>  
> and get_first_[line/quad/hex]_index() 
> <https://dealii.org/developer/doxygen/deal.II/classFiniteElementData.html#a078fb199880f2b29fc72d9c01d317bd0>
>  
> functions, which does suggest that I'm wrong about this, so I'd be happy to 
> be corrected on this point.  (Question: What does a vertex/edge/face DoF 
> even mean for a DG element?) Maybe these are the exact functions that you 
> were looking for (and maybe with this , if the [vertex,edge,face,cell] 
> groupings do apply, you can map from the cell DoF->associated geometric 
> entity using these functions). But I don't that they would help 
> understanding the ordering (e.g. lexiographic or something else; element 
> base index fast or spacedim fast, etc.) within each grouping -- again, 
> that's probably a decision that was left to the designer of each FE class, 
> and I don't think that there's been a need for any consistency between the 
> classes. But again, I'm just putting my thoughts out there and anticipate 
> being corrected on these points.
>

In my humble opinion the FE should only be implemented in terms of a 
reference element. This is, as far as I can tell, the case for all elements 
in Deal.II. The mapping though is problematic on complicated domains and I 
saw that one needs to renumber or switch signs for some shape functions if 
one wants to satisfy continuity conditions for the respective spaces. This 
is partly taken care of by the DoFHandler (see renumbering of FE_Q on 
faces) class which confused me, while other parts like the sign change is 
not. I guess I got the idea. Also, the concept of (generalized) support 
points is nice but sometimes it is not really clear to me if this is 
helpful in the context of element pairings (for example in the language of 
finite element exterior calculus where this concept does not really make 
sense globally). I may be wrong and miss a point here  and would love to 
get a lesson.
 

>
> I get the feeling that there is no overlap in assignment of DoFs between 
> geometric entities. Consider each of these entities not to have closure: so 
> you're actually referring to DoFs associated with the cell interior (excl. 
> faces, edges and vertices), face interior (excl edges and vertices) and 
> edge interiors (excl. vertices). To me, that the only sensible 
> interpretation of this. 
>
 
This is partly what is not fully clear to me. Vertex/edge/face/volume dofs 
(again FEEC) depends on the functionals that define them and less on 
location (a Raviart-Thomas or Nedelec shape function is not regular enough 
for point evaluations, although viewing them as a polynomial suggests so). 


> Thanks for giving this so much thought, and for the interesting questions. 
> I've certainly go

[deal.II] Re: Mean Value Constraints

[Konrad Simon]

Little correction: I wrote "In the vmult() function remove the mean value 
(i.e., project the rhs on the orthogonal complement of the kernel of the 
kernel)", I meant "In the vmult() function remove the mean value after you 
multiply (i.e., project the rhs on the orthogonal complement of the kernel 
of the system matrix)"

This paper <https://epubs.siam.org/doi/abs/10.1137/S0036144503426074> by 
Bochev and Lehoucq describes the technique with a toy model (Poisson with 
purely natural BCs). It was not new at that time but to the best of my 
knowledge they were the first ones to write this down in a proper way. The 
technique I mentioned above works for other kernels in system matrices as 
well. For example in linearized elasticity one maybe needs to remove not 
the mean value but rigid motions (well not really rigid but infinitesimal 
rotations - invariance under rotations cannot be linear) from the kernel. 
Such motions do not affect stress but constitute of nontrivial 
displacements. Also in this case one can project the rhs onto the 
orthogonal complement of the kernel of the system matrix. In this case the 
kernel has dimension 3 in 2D (one for scaling+90 degree rotation and two 
independent shifts) and dimension 6 in 3D. Once you project the rhs a 
Krylov solver can deal with your singular problem. 

Cheers,
Konrad 
On Saturday, December 26, 2020 at 11:06:36 AM UTC+1 Konrad Simon wrote:

> Hi,
>
> On Saturday, December 26, 2020 at 6:43:21 AM UTC+1 smetca...@gmail.com 
> wrote:
>
>> Hi all,
>>
>> Does anyone here have any experience applying mean value constraints 
>> (specifically with periodic boundary conditions)? I'm having some trouble. 
>> As far as I can tell, there are two approaches (both of which I have been 
>> able to get working). The first approach is to just add the mean value 
>> constraint straight into the constraint matrix. I got this working just 
>> fine and was able to get an output but the problem is that this condition 
>> seems to make my matrices dense(?) or, at the very least, require a huge 
>> amount of memory hindering this approach for the types of fine spatial 
>> meshes I require for this application. 
>>
>
> If you think of the mean value constraint as an additional row it looks 
> like (for Poisson equation) u_0+u_1+...+u_n=0 and hence it couples all 
> DoFs. This is mathematically ok but not efficient since this constraint is 
> non-local. What one can also do is just constrain one DoF to a specific 
> value (this would also remove rigid motion in elasticity). But think about 
> your solution variable: If it is in the Sobolev space H^1 then point 
> evaluations may not be defined for dimension larger than 2. Similarly if, 
> for example, the pressure in a mechanical or fluid problem is often just in 
> L^2. Point evaluations do not make sense there at all. 
>  
>
>> The other approach is to just apply periodic boundary conditions to the 
>> matrices, try to solve the system anyway and, because UMFPACK is awesome, 
>> somehow manage to get a solution out of all of this (which will of course 
>> only be defined up to a constant). The problem here is that the solution I 
>> am getting is of order 10^12! When I postprocess and subtract the mean 
>> value, what I get looks correct but I am losing a lot of accuracy in the 
>> process but it does work for the fine meshes that I require. Anyone have 
>> any ideas or do I just have to live with the reduced accuracy here? Thanks!
>>
>
> This is of course like solving an ill-posed problem but some direct 
> solvers are smart enough to deal with zero pivots (redundant rows).
>
> Both approaches above do work but only for small and non-parallel 
> problems. This is what you can do:
>
> Many people do not think of the fact that even iterative solvers can deal 
> with singular problems if the right hand side is in the orthogonal 
> complement of the kernel. For the mean value constraint do the following:
>
> 1.  Setup boundary conditions in constraints if necessary. 
> 2. Assemble the system matrix and rhs and apply your constraints
> 3. Create a class that wraps your system matrix. This class must provide a 
> vmult() function
> 4. Before starting the iterative solver remove the mean value, see step-57 
> <https://www.dealii.org/current/doxygen/deal.II/step_56.html#StokesProblemprocess_solution>
>  
> 4. In the vmult() function remove the mean value (i.e., project the rhs on 
> the orthogonal complement of the kernel of the kernel)
> 5. use the wrapped matrix in the solver.
>
> 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

[deal.II] Re: Mean Value Constraints

[Konrad Simon]

Hi,

On Saturday, December 26, 2020 at 6:43:21 AM UTC+1 smetca...@gmail.com 
wrote:

> Hi all,
>
> Does anyone here have any experience applying mean value constraints 
> (specifically with periodic boundary conditions)? I'm having some trouble. 
> As far as I can tell, there are two approaches (both of which I have been 
> able to get working). The first approach is to just add the mean value 
> constraint straight into the constraint matrix. I got this working just 
> fine and was able to get an output but the problem is that this condition 
> seems to make my matrices dense(?) or, at the very least, require a huge 
> amount of memory hindering this approach for the types of fine spatial 
> meshes I require for this application. 
>

If you think of the mean value constraint as an additional row it looks 
like (for Poisson equation) u_0+u_1+...+u_n=0 and hence it couples all 
DoFs. This is mathematically ok but not efficient since this constraint is 
non-local. What one can also do is just constrain one DoF to a specific 
value (this would also remove rigid motion in elasticity). But think about 
your solution variable: If it is in the Sobolev space H^1 then point 
evaluations may not be defined for dimension larger than 2. Similarly if, 
for example, the pressure in a mechanical or fluid problem is often just in 
L^2. Point evaluations do not make sense there at all. 
 

> The other approach is to just apply periodic boundary conditions to the 
> matrices, try to solve the system anyway and, because UMFPACK is awesome, 
> somehow manage to get a solution out of all of this (which will of course 
> only be defined up to a constant). The problem here is that the solution I 
> am getting is of order 10^12! When I postprocess and subtract the mean 
> value, what I get looks correct but I am losing a lot of accuracy in the 
> process but it does work for the fine meshes that I require. Anyone have 
> any ideas or do I just have to live with the reduced accuracy here? Thanks!
>

This is of course like solving an ill-posed problem but some direct solvers 
are smart enough to deal with zero pivots (redundant rows).

Both approaches above do work but only for small and non-parallel problems. 
This is what you can do:

Many people do not think of the fact that even iterative solvers can deal 
with singular problems if the right hand side is in the orthogonal 
complement of the kernel. For the mean value constraint do the following:

1.  Setup boundary conditions in constraints if necessary. 
2. Assemble the system matrix and rhs and apply your constraints
3. Create a class that wraps your system matrix. This class must provide a 
vmult() function
4. Before starting the iterative solver remove the mean value, see step-57 
<https://www.dealii.org/current/doxygen/deal.II/step_56.html#StokesProblemprocess_solution>
 
4. In the vmult() function remove the mean value (i.e., project the rhs on 
the orthogonal complement of the kernel of the kernel)
5. use the wrapped matrix in the solver.

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/67679343-a4c9-46c3-96b9-dd037445046fn%40googlegroups.com.

Re: [deal.II] Div-conforming elements in 3D

[Konrad Simon]

Dear Deal.II community,

Suppose I chose FE_Nedelec<3> of some order as finite element.

1. Is there a way to query for a given cell_dof_index whether 
it is a face_dof and to get the face_index?

2. If this dof is a line_dof is it by definition also a 
face_dof? I need to exclude line_dofs and treat them 
separately.

3. Is there a convention that (if I am not using FESystem) in the set of 
all cell_dofs all line_dofs come before all 
face_dofs and only then the volume_dofs?

I am asking since I am trying to come up with a patch for RaviartThomas and 
a pattern to implement it in the current structures without interfering 
with them (too much).

Best,
Konrad
On Wednesday, December 16, 2020 at 9:33:57 PM UTC+1 Konrad Simon wrote:

> Dear Jean-Paul, dear Deal.II community,
>
> I partially solved the problem of sign flipping and permuting the degrees 
> of freedom and would like to come up with a merge-request at least for the 
> RaviartThomas elements. I have some questions before I can go on and 
> guidance would be appreciated.
>
> I decided to go a similar way that is taken for FE_Q elements. There one 
> has only a permutation of dofs on non-standard or flipped or rotated faces 
> since a sign change in the orientation does not matter. A vector of size 
> n_faces_per_cell of which each contains a table mapping each face_dof and 
> each of the 8 combinations of the flags (non-standard | flipped | rotated) 
> to an offset of that face_dof. The system behind that relies on the 
> lexicographic ordering of dofs if I understand correctly. My questions:
>
> 1. Is there a similar numbering scheme for dofs on faces for RT elements. 
> There one has face moments as dofs and it is not fully clear to me what 
> lexicographic means there (if such a principle is usable there at all).
>
> 2. If the order (and sign) of face_dofs if permuted for each of the flags 
> (non-standard | flipped | rotated) then in what order are the permutations 
> applied? (They do not commute)
>
> I would like to avoid hard-coding that (I think up to order 2 would be 
> doable -> but already 54 face_dofs). 
>
> I believe the table similar 
> to adjust_quad_dof_index_for_face_orientation_table of FE_Q_Base makes 
> sense but I guess every FE inheriting from FE_PolyTensor must implement 
> such a table then. This is a lot of work and more complicated for Nedelec 
> elements.
>
> Best,
> Konrad
> On Saturday, December 12, 2020 at 8:07:19 PM UTC+1 Konrad Simon wrote:
>
>> Hi Jean-Paul,
>>
>> On Thursday, December 10, 2020 at 11:39:08 PM UTC+1 Jean-Paul Pelteret 
>> wrote:
>>
>>> HI Konrad,
>>>
>>> I have no familiarity with the H-div elements, so I could be wrong with 
>>> this suggestion...
>>>
>>> The Fe_Nedelec element suffered from a similar issue, where adjacent 
>>> cells had to agree on the edge sense/directionality. This requirement could 
>>> not be unconditionally guaranteed for that element type, hence the 
>>> following note in its introduction 
>>> <https://dealii.org/current/doxygen/deal.II/classFE__Nedelec.html>
>>>
>>> Even if this element is implemented for two and three space dimensions, 
>>> the definition of the node values relies on consistently oriented faces in 
>>> 3D. Therefore, care should be taken on complicated meshes.
>>>
>>
>> Yes, this issue is similar. I assume I will need to check that for these 
>> elements as well since I also need them.
>>
>>
>>> The more recent FE_NedelecSZ 
>>> <https://dealii.org/current/doxygen/deal.II/classFE__NedelecSZ.html> 
>>> element resolves the issue by querying the vertex numbers that make up each 
>>> edge/face and applying a rule that gives each face and its bounding edges a 
>>> direction. This means that the orientation of the face/edges is always 
>>> consistent no matter which of the two cells that share the face you're on. 
>>> The details are more explicitly laid out in the FE_NedelecSZ 
>>> <https://dealii.org/current/doxygen/deal.II/classFE__NedelecSZ.html> 
>>> introduction, as well as the two papers that are listed therein. 
>>>
>>
>> This is a nice alternative for the Nedelec element. Often though such 
>> elements are used in conjunction with H1 or H(div) conformal elements and 
>> need to satisfy a compatibility condition (stronger than just LBB, i.e., 
>> they must form a bounded co-chain complex that is a subcomplex to something 
>> else). I am not sure that this element does satisfy this condition together 
>> with classical RT elements etc.
>>
>>
>>> Now, I'm by no means certain tha

Re: [deal.II] Div-conforming elements in 3D

[Konrad Simon]

Dear Jean-Paul, dear Deal.II community,

I partially solved the problem of sign flipping and permuting the degrees 
of freedom and would like to come up with a merge-request at least for the 
RaviartThomas elements. I have some questions before I can go on and 
guidance would be appreciated.

I decided to go a similar way that is taken for FE_Q elements. There one 
has only a permutation of dofs on non-standard or flipped or rotated faces 
since a sign change in the orientation does not matter. A vector of size 
n_faces_per_cell of which each contains a table mapping each face_dof and 
each of the 8 combinations of the flags (non-standard | flipped | rotated) 
to an offset of that face_dof. The system behind that relies on the 
lexicographic ordering of dofs if I understand correctly. My questions:

1. Is there a similar numbering scheme for dofs on faces for RT elements. 
There one has face moments as dofs and it is not fully clear to me what 
lexicographic means there (if such a principle is usable there at all).

2. If the order (and sign) of face_dofs if permuted for each of the flags 
(non-standard | flipped | rotated) then in what order are the permutations 
applied? (They do not commute)

I would like to avoid hard-coding that (I think up to order 2 would be 
doable -> but already 54 face_dofs). 

I believe the table similar 
to adjust_quad_dof_index_for_face_orientation_table of FE_Q_Base makes 
sense but I guess every FE inheriting from FE_PolyTensor must implement 
such a table then. This is a lot of work and more complicated for Nedelec 
elements.

Best,
Konrad
On Saturday, December 12, 2020 at 8:07:19 PM UTC+1 Konrad Simon wrote:

> Hi Jean-Paul,
>
> On Thursday, December 10, 2020 at 11:39:08 PM UTC+1 Jean-Paul Pelteret 
> wrote:
>
>> HI Konrad,
>>
>> I have no familiarity with the H-div elements, so I could be wrong with 
>> this suggestion...
>>
>> The Fe_Nedelec element suffered from a similar issue, where adjacent 
>> cells had to agree on the edge sense/directionality. This requirement could 
>> not be unconditionally guaranteed for that element type, hence the 
>> following note in its introduction 
>> <https://dealii.org/current/doxygen/deal.II/classFE__Nedelec.html>
>>
>> Even if this element is implemented for two and three space dimensions, 
>> the definition of the node values relies on consistently oriented faces in 
>> 3D. Therefore, care should be taken on complicated meshes.
>>
>
> Yes, this issue is similar. I assume I will need to check that for these 
> elements as well since I also need them.
>
>
>> The more recent FE_NedelecSZ 
>> <https://dealii.org/current/doxygen/deal.II/classFE__NedelecSZ.html> 
>> element resolves the issue by querying the vertex numbers that make up each 
>> edge/face and applying a rule that gives each face and its bounding edges a 
>> direction. This means that the orientation of the face/edges is always 
>> consistent no matter which of the two cells that share the face you're on. 
>> The details are more explicitly laid out in the FE_NedelecSZ 
>> <https://dealii.org/current/doxygen/deal.II/classFE__NedelecSZ.html> 
>> introduction, as well as the two papers that are listed therein. 
>>
>
> This is a nice alternative for the Nedelec element. Often though such 
> elements are used in conjunction with H1 or H(div) conformal elements and 
> need to satisfy a compatibility condition (stronger than just LBB, i.e., 
> they must form a bounded co-chain complex that is a subcomplex to something 
> else). I am not sure that this element does satisfy this condition together 
> with classical RT elements etc.
>
>
>> Now, I'm by no means certain that underlying issue that you're seeing 
>> here is the same as that which plagued the canonical Nedelec element 
>> (although your comment #3 makes me suspect that it could be). But if it is, 
>> then perhaps what was done to fix the orientation conflict in the Zaglmayr 
>> formulation of the Nedelec element can serve as inspiration for a fix to 
>> the BDM and RT elements?
>>
>> Thanks for looking into this, and I hope that you have some success at 
>> finding a solution to the problem. Naturally, if there's anything that we 
>> can do to help please let us know. It would be very nice to have a robust 
>> implementation to these elements in the library, so anything that you might 
>> be willing to contribute would be greatly appreciated!
>>
>  
> Lowest order RaviartThomas is already fixed. I will hunt down this issue 
> as good as I can. I might need some help with local and global dof 
> numbering conventions. I will also try to move the discussion to the issue 
> page <https://github.com/dea

Re: [deal.II] Div-conforming elements in 3D

[Konrad Simon]

Hi Jean-Paul,

On Thursday, December 10, 2020 at 11:39:08 PM UTC+1 Jean-Paul Pelteret 
wrote:

> HI Konrad,
>
> I have no familiarity with the H-div elements, so I could be wrong with 
> this suggestion...
>
> The Fe_Nedelec element suffered from a similar issue, where adjacent cells 
> had to agree on the edge sense/directionality. This requirement could not 
> be unconditionally guaranteed for that element type, hence the following note 
> in its introduction 
> <https://dealii.org/current/doxygen/deal.II/classFE__Nedelec.html>
>
> Even if this element is implemented for two and three space dimensions, 
> the definition of the node values relies on consistently oriented faces in 
> 3D. Therefore, care should be taken on complicated meshes.
>

Yes, this issue is similar. I assume I will need to check that for these 
elements as well since I also need them.


> The more recent FE_NedelecSZ 
> <https://dealii.org/current/doxygen/deal.II/classFE__NedelecSZ.html> 
> element resolves the issue by querying the vertex numbers that make up each 
> edge/face and applying a rule that gives each face and its bounding edges a 
> direction. This means that the orientation of the face/edges is always 
> consistent no matter which of the two cells that share the face you're on. 
> The details are more explicitly laid out in the FE_NedelecSZ 
> <https://dealii.org/current/doxygen/deal.II/classFE__NedelecSZ.html> 
> introduction, as well as the two papers that are listed therein. 
>

This is a nice alternative for the Nedelec element. Often though such 
elements are used in conjunction with H1 or H(div) conformal elements and 
need to satisfy a compatibility condition (stronger than just LBB, i.e., 
they must form a bounded co-chain complex that is a subcomplex to something 
else). I am not sure that this element does satisfy this condition together 
with classical RT elements etc.


> Now, I'm by no means certain that underlying issue that you're seeing here 
> is the same as that which plagued the canonical Nedelec element (although 
> your comment #3 makes me suspect that it could be). But if it is, then 
> perhaps what was done to fix the orientation conflict in the Zaglmayr 
> formulation of the Nedelec element can serve as inspiration for a fix to 
> the BDM and RT elements?
>
> Thanks for looking into this, and I hope that you have some success at 
> finding a solution to the problem. Naturally, if there's anything that we 
> can do to help please let us know. It would be very nice to have a robust 
> implementation to these elements in the library, so anything that you might 
> be willing to contribute would be greatly appreciated!
>
 
Lowest order RaviartThomas is already fixed. I will hunt down this issue as 
good as I can. I might need some help with local and global dof numbering 
conventions. I will also try to move the discussion to the issue page 
<https://github.com/dealii/dealii/issues/7970> on github.

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/90968d76-09ae-4dd5-82a8-c49bc6383792n%40googlegroups.com.

Re: [deal.II] Div-conforming elements in 3D

[Konrad Simon]

Dear David, dear all,

It seems like the issue is much more involved than I initially thought. 
Some updates since I started going into this issue 
<https://github.com/dealii/dealii/issues/7970>.

1. I could solve the problem for lowest order RaviartThomas Elements (RT). 
This is already good.

2. The fix does not work for lowest order BDM elements.

3. The fix does not work for higher order elements, neither BDM nor RT 
elements
Suspected reason: For shape functions of higher order elements (whose 
normal components) are supported on a face flipping the sign on a flipped 
face is not enough since the flipped face changes/reverts the ordering of 
the face dofs hence the translation of local to global numbering of degrees 
of freedom changes for cells with flipped faces. I am not sure but at least 
I suspect that. Is this being taken care of by Deal.II somewhere?

I will go on chasing the bug and hopefully come up with a merge request 
soon.

Best,
Konrad

On Wednesday, December 9, 2020 at 3:44:58 PM UTC+1 Konrad Simon wrote:

> Hi David, 
>
> Many thanks for the hint. After some research I believe I stumbled over this 
> issue#7970 <https://github.com/dealii/dealii/issues/7970>? Let's see what 
> I can do.
>
> Best,
> Konrad
> On Tuesday, December 8, 2020 at 10:56:44 PM UTC+1 Wells, David wrote:
>
>> Hi Konrad,
>>
>> This isn't a very satisfying answer but, to the best of my knowledge, 
>> these Hdiv elements weren't written with general unstructured 3D meshes in 
>> mind. I suspect, without looking more closely, that this is a problem on 
>> our side.
>>
>> We document this for the FE_RaviartThomasNodal class but for 
>> FE_RaviartThomas this is just a comment in the source file. Would you be 
>> willing to write a patch explaining things better for the other Hdiv 
>> classes?
>>
>> Best,
>> David
>> --
>> *From:* dea...@googlegroups.com  on behalf of 
>> Konrad Simon 
>> *Sent:* Monday, December 7, 2020 12:48 PM
>> *To:* deal.II User Group 
>> *Subject:* [deal.II] Div-conforming elements in 3D 
>>  
>> Dear Deal.ii community, 
>>
>> Sorry for spamming the mailing list again. I posted some of the below 
>> stuff already in another post but that might have been with too many 
>> unnecessary details and confusing ideas from my side.
>>
>> I am not sure if I stumbled over something and would like to report on 
>> that. 
>>
>> Goal is to project a given function onto a finite element space over a 
>> certain mesh. The focus here shall be on H(div)-conforming elements such as 
>> Raviart-Thomas and Brezzi-Douglas-Marini Elements.
>>
>> I do so with several different meshes and I noticed an inconsistency in 
>> the projection on meshes with, say, not so simple topology (e.g., holes and 
>> voids). 
>>
>> Observations:
>> 1. The problematic meshes do have faces with non-standard orientation. My 
>> first thought was this was not carried over to the mapping of the 
>> H(div)-basis but I was wrong. The mesh orientation is consistent (as far as 
>> I can tell). A cube, for example, is not problematic. Please see two 
>> example plots below.
>>
>> 2. The issue arises with both Raviart-Thomas and Brezzi-Douglas-Marini 
>> Elements.
>>
>> 3. The issue arises in 3D only. 
>>
>> Any ideas what I ran into?
>>
>> Best,
>> Konrad[image: plate_with_hole.png][image: hyper_shell.png]
>>
>> -- 
>> 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+un...@googlegroups.com.
>> To view this discussion on the web visit 
>> https://groups.google.com/d/msgid/dealii/fa304b41-f83f-4c96-a2f7-c4e9e813553fn%40googlegroups.com
>>  
>> <https://groups.google.com/d/msgid/dealii/fa304b41-f83f-4c96-a2f7-c4e9e813553fn%40googlegroups.com?utm_medium=email_source=footer>
>> .
>>
>

-- 
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/883fee09-6fba-4d38-914d-869e2e7058adn%40googlegroups.com.

Re: [deal.II] Div-conforming elements in 3D

[Konrad Simon]

Hi David, 

Many thanks for the hint. After some research I believe I stumbled over this 
issue#7970 <https://github.com/dealii/dealii/issues/7970>? Let's see what I 
can do.

Best,
Konrad
On Tuesday, December 8, 2020 at 10:56:44 PM UTC+1 Wells, David wrote:

> Hi Konrad,
>
> This isn't a very satisfying answer but, to the best of my knowledge, 
> these Hdiv elements weren't written with general unstructured 3D meshes in 
> mind. I suspect, without looking more closely, that this is a problem on 
> our side.
>
> We document this for the FE_RaviartThomasNodal class but for 
> FE_RaviartThomas this is just a comment in the source file. Would you be 
> willing to write a patch explaining things better for the other Hdiv 
> classes?
>
> Best,
> David
> --
> *From:* dea...@googlegroups.com  on behalf of 
> Konrad Simon 
> *Sent:* Monday, December 7, 2020 12:48 PM
> *To:* deal.II User Group 
> *Subject:* [deal.II] Div-conforming elements in 3D 
>  
> Dear Deal.ii community, 
>
> Sorry for spamming the mailing list again. I posted some of the below 
> stuff already in another post but that might have been with too many 
> unnecessary details and confusing ideas from my side.
>
> I am not sure if I stumbled over something and would like to report on 
> that. 
>
> Goal is to project a given function onto a finite element space over a 
> certain mesh. The focus here shall be on H(div)-conforming elements such as 
> Raviart-Thomas and Brezzi-Douglas-Marini Elements.
>
> I do so with several different meshes and I noticed an inconsistency in 
> the projection on meshes with, say, not so simple topology (e.g., holes and 
> voids). 
>
> Observations:
> 1. The problematic meshes do have faces with non-standard orientation. My 
> first thought was this was not carried over to the mapping of the 
> H(div)-basis but I was wrong. The mesh orientation is consistent (as far as 
> I can tell). A cube, for example, is not problematic. Please see two 
> example plots below.
>
> 2. The issue arises with both Raviart-Thomas and Brezzi-Douglas-Marini 
> Elements.
>
> 3. The issue arises in 3D only. 
>
> Any ideas what I ran into?
>
> Best,
> Konrad[image: plate_with_hole.png][image: hyper_shell.png]
>
> -- 
> 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+un...@googlegroups.com.
> To view this discussion on the web visit 
> https://groups.google.com/d/msgid/dealii/fa304b41-f83f-4c96-a2f7-c4e9e813553fn%40googlegroups.com
>  
> <https://groups.google.com/d/msgid/dealii/fa304b41-f83f-4c96-a2f7-c4e9e813553fn%40googlegroups.com?utm_medium=email_source=footer>
> .
>

-- 
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/2570e306-91b1-4d03-90fa-9a3777066264n%40googlegroups.com.

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 
<https://dealii.org/current/doxygen/deal.II/DEALGlossary.html#GlossFaceOrientation>
 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.

[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(n

Re: [deal.II] Tools for parallel debugging

[Konrad Simon]

Thank you, Wolfgang and Daniel. Seems like I will go with the command line 
then. I was just wondering if people here use Eclipse's PTP which sounded 
like a good graphical tool. In my case it frequently crashes or simply gets 
stuck.

Best,
Konrad

On Sunday, October 11, 2020 at 11:35:34 PM UTC+2 Wolfgang Bangerth wrote:

> On 10/11/20 3:26 PM, Daniel Arndt wrote:
> > 
> > have a look at 
> > 
> https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#how-do-i-debug-mpi-programs
>  
> > <
> https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fdealii%2Fdealii%2Fwiki%2FFrequently-Asked-Questions%23how-do-i-debug-mpi-programs=02%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cab82fabba6ff45f8e49708d86e2c6359%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637380484218698974=jECVtixHfq%2BTSXkb1xww9SYWSyGwOonjWg%2Fqug7zLPk%3D=0>.
>  
>
> > For me, the "mpiexec -np n gdb ./executable" approach normally works 
> fairly well.
> > 
>
> There is also an interactive demonstration of this in lecture 41.25:
> https://www.math.colostate.edu/~bangerth/videos.676.41.25.html
>
> Best
> W.
>
> -- 
> 
> Wolfgang Bangerth email: bang...@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/3ed1b810-e78e-4cf0-9712-70aec571875en%40googlegroups.com.

Re: [deal.II] Tools for parallel debugging

[Konrad Simon]

Oh, seems like I missed that. Thanks you!

Konrad

On Sunday, October 11, 2020 at 11:35:34 PM UTC+2 Wolfgang Bangerth wrote:

> On 10/11/20 3:26 PM, Daniel Arndt wrote:
> > 
> > have a look at 
> > 
> https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#how-do-i-debug-mpi-programs
>  
> > <
> https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fdealii%2Fdealii%2Fwiki%2FFrequently-Asked-Questions%23how-do-i-debug-mpi-programs=02%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cab82fabba6ff45f8e49708d86e2c6359%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637380484218698974=jECVtixHfq%2BTSXkb1xww9SYWSyGwOonjWg%2Fqug7zLPk%3D=0>.
>  
>
> > For me, the "mpiexec -np n gdb ./executable" approach normally works 
> fairly well.
> > 
>
> There is also an interactive demonstration of this in lecture 41.25:
> https://www.math.colostate.edu/~bangerth/videos.676.41.25.html
>
> Best
> W.
>
> -- 
> 
> Wolfgang Bangerth email: bang...@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/3af7489d-9b56-4796-af4c-498291829c2an%40googlegroups.com.

[deal.II] Tools for parallel debugging

[Konrad Simon]

Hi deal.ii community,

I have a short question about tools (which probably also concerns many 
other people). 

I am using eclipse for C++ development and I frequently use the debugger. 
So far everything fine. Now I need to use parallel debugging tools for MPI 
but my eclipse crashes many times and/or is super slow.

What tools are you using? Any recommendations or hints?

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/d89c4b39-a069-4825-b96b-741665e55095n%40googlegroups.com.

Re: [deal.II] Re: Evaluating FE-solution on distributed mesh, semi-Lagrangian method

[Konrad Simon]

Dear all,

I finally had time to go back to my code. Quick problem summary: I am 
currently working on a semi-Lagrangian method for an advection-diffusion 
equation with a multiscale method that I designed for advection dominated 
processes. During runtime I must evaluate the (known) multiscale 
FE-solution object at a previous time step. 

Problem: The mesh is distributed and the object that will give me the 
evaluated data at a point p is possibly owned on another processor. The 
processor that owes the object is the same as the one owing the cell that 
surrounds p. To find the correct cell on the same processor that needs the 
information I am doing things similarly as the FEFieldFunction class to 
find the cell. But what if p is in a cell that is not locally owned? It 
could be in a ghost cell or an artificial one (also ghosts no not owe the 
data I need).

Now, I have have to manually MPI communicate such points to find out the 
MPI rank that owes the cell surrounding p as I found out with your help. 
Any way of doing this efficiently? I would like to avoid the brute force 
version: communicating all points to all ranks and search the owned cells. 
(Hint: Mostly (but not always) p will be in a ghost layer. )

Any way to get the actual MPI rank of the owner?

Any help would be much appreciated.

Best,
Konrad

On Thursday, May 14, 2020 at 4:10:33 PM UTC+2 Konrad Simon wrote:

> Thank you, Bruno. :-)
>
>
> On Tuesday, May 12, 2020 at 2:50:35 PM UTC+2, Bruno Turcksin wrote:
>
>> Konrad, 
>>
>> There is nothing out of the box. However, deal.II uses p4est which can 
>> use more that one layer of ghost cells. So you should take a look 
>> there to see how hard it is to change the code. If that's the way you 
>> want to go, I am sure that someone will be able to provide you with 
>> some guidance. 
>>
>> Best, 
>>
>> Bruno 
>>
>> Le mar. 12 mai 2020 à 06:04, Konrad Simon  a écrit : 
>> > 
>> > Thank you, Bruno. What is sometimes done for semi-Lagrangian methods is 
>> that one defines a halo region around the locally owned cells (like ghost 
>> cells) that contains information from previous time steps and then limits 
>> the time step so that one never leaves the halo region when asking for 
>> information from previous time steps. Is there a way in Deal.ii to control 
>> the size of the halo region. In my case I have pretty uniform grids if it 
>> helps. 
>> > 
>> > Best, 
>> > Konrad 
>> > 
>> > On Monday, May 11, 2020 at 2:49:21 PM UTC+2, Bruno Turcksin wrote: 
>> >> 
>> >> Konrad, 
>> >> 
>> >> Unfortunately, you will need to do the communication yourself. You can 
>> only evaluate the solution on cells that locally owned or on ghost cells. 
>> >> 
>> >> Best, 
>> >> 
>> >> Bruno 
>> > 
>> > -- 
>> > 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 a topic in the 
>> Google Groups "deal.II User Group" group. 
>> > To unsubscribe from this topic, visit 
>> https://groups.google.com/d/topic/dealii/9li9twVxrGE/unsubscribe. 
>>
> > To unsubscribe from this group and all its topics, send an email to 
>> dea...@googlegroups.com. 
>>
> > To view this discussion on the web visit 
>> https://groups.google.com/d/msgid/dealii/bfc6dd81-abb3-421c-8f6c-539c234b624f%40googlegroups.com.
>>  
>>
>>
>

-- 
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/1b60b7a8-361a-4f2f-9216-50190ab7e950n%40googlegroups.com.

Re: [deal.II] Re: Evaluating FE-solution on distributed mesh, semi-Lagrangian method

[Konrad Simon]

Thank you, Bruno. :-)

On Tuesday, May 12, 2020 at 2:50:35 PM UTC+2, Bruno Turcksin wrote:
>
> Konrad, 
>
> There is nothing out of the box. However, deal.II uses p4est which can 
> use more that one layer of ghost cells. So you should take a look 
> there to see how hard it is to change the code. If that's the way you 
> want to go, I am sure that someone will be able to provide you with 
> some guidance. 
>
> Best, 
>
> Bruno 
>
> Le mar. 12 mai 2020 à 06:04, Konrad Simon  > a écrit : 
> > 
> > Thank you, Bruno. What is sometimes done for semi-Lagrangian methods is 
> that one defines a halo region around the locally owned cells (like ghost 
> cells) that contains information from previous time steps and then limits 
> the time step so that one never leaves the halo region when asking for 
> information from previous time steps. Is there a way in Deal.ii to control 
> the size of the halo region. In my case I have pretty uniform grids if it 
> helps. 
> > 
> > Best, 
> > Konrad 
> > 
> > On Monday, May 11, 2020 at 2:49:21 PM UTC+2, Bruno Turcksin wrote: 
> >> 
> >> Konrad, 
> >> 
> >> Unfortunately, you will need to do the communication yourself. You can 
> only evaluate the solution on cells that locally owned or on ghost cells. 
> >> 
> >> Best, 
> >> 
> >> Bruno 
> > 
> > -- 
> > 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 a topic in the 
> Google Groups "deal.II User Group" group. 
> > To unsubscribe from this topic, visit 
> https://groups.google.com/d/topic/dealii/9li9twVxrGE/unsubscribe. 
> > To unsubscribe from this group and all its topics, send an email to 
> dea...@googlegroups.com . 
> > To view this discussion on the web visit 
> https://groups.google.com/d/msgid/dealii/bfc6dd81-abb3-421c-8f6c-539c234b624f%40googlegroups.com.
>  
>
>

-- 
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/f7253b78-0633-4a35-abf2-ea8c76576d71%40googlegroups.com.

[deal.II] Evaluating FE-solution on distributed mesh, semi-Lagrangian method

[Konrad Simon]

Dear all,

I am currently working on a semi-Lagrangian method for an 
advection-diffusion equation. During runtime I must evaluate the (known) 
FE-solution at a previous time step but the mesh is distributed. 

Problem: It can well be that I must evaluate the solution at a point that 
is in a cell that is not locally owed by the processor requesting that 
value or not even a ghost cell. I am aware of the FEFieldFunction class and 
that it says that an exception will be thrown if the point is found to be 
in an artificial cell. 

Now, does anyone here have experience with that? Am I doomed and will I 
have to manually MPI communicate such points to any processor and test if 
they own the relevant cell?

Any help would be much appreciated.

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/9e19c2b8-29e2-48fd-bc99-18b912696cfd%40googlegroups.com.

Re: [deal.II] Installation with spack - no access to adol-c

[Konrad Simon]

Jean-Paul, thanks a lot. Works! :-)

Cheers,
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/47803471-9b2b-4af4-aa9c-b77d39a68af6%40googlegroups.com.

[deal.II] Installation with spack - no access to adol-c

[Konrad Simon]

Hi all,

A colleague and myself had problems to install dealii 9.1.1 via spack.

Quick question: Did anyone have problems to install deal.ii via spack 
lately? I see that it depends on the development version of adol-c but this 
one requires a gitlab account and access to adol-c. Does deal.ii v9.1.1 
really depend on the development version of adol-c? Other versions <= 2.6.3 
are on github.

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/a09a369b-ef5c-48d4-bdf5-1c3c681d614c%40googlegroups.com.

[deal.II] Re: Interpolating to globally refined distributed mesh

[Konrad Simon]

Dear deal.ii community,

I solved it, sorry for bugging you with it but simple mistakes can bug you 
for long... 

I simply forgot to re-distribute the dofs for the finite element after 
refining the mesh. Ooofff :-/

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/00bb0838-ea8a-4704-8a13-94089509a3c6%40googlegroups.com.

[deal.II] Re: Interpolating to globally refined distributed mesh

[Konrad Simon]

Add-on: The problem is in the function 
RefineInterpolate::refine_and_interpolate_on_distributed_mesh()

-- 
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/51923469-509d-4c73-ac27-cc5e5b79cce1%40googlegroups.com.

[deal.II] Re: Interpolating to globally refined distributed mesh

[Konrad Simon]

Add-on: The problem is in the function 
RefineInterpolate::refine_and_interpolate_on_distributed_mesh()

-- 
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/beca16b0-92ec-46f3-90bf-a1ccfc2913fb%40googlegroups.com.

[deal.II] Interpolating to globally refined distributed mesh

[Konrad Simon]

Dear deal.ii community,

I am trying to interpolate a distributed solution vector onto a globally 
refined distributed mesh. Actually I need to globally refine more than once 
but I believe I can iterate this procedure (is there a more efficient way? 
The finite element can be vector or scalar valued, and the solution vector 
a block vector)

However, despite looking at tutorials and the documentation I get an error 
(segmentation fault) without further explanation. I did not manage to track 
down the error - probably I am using the interface not correctly. Anyone 
has experience with this?

A minimal example is attached. (I used deal.ii 9.2.0 but it should work 
with 9.1.1, too.)

Best regards,
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/9978ebc5-4ff9-411f-b892-d8220a2146ec%40googlegroups.com.
##
#  CMake script for the step-40 tutorial program:
##

# Set the name of the project and target:
SET(TARGET "interpolate")

# Declare all source files the target consists of. Here, this is only
# the one step-X.cc file, but as you expand your project you may wish
# to add other source files as well. If your project becomes much larger,
# you may want to either replace the following statement by something like
#  FILE(GLOB_RECURSE TARGET_SRC  "source/*.cc")
#  FILE(GLOB_RECURSE TARGET_INC  "include/*.h")
#  SET(TARGET_SRC ${TARGET_SRC}  ${TARGET_INC})
# or switch altogether to the large project CMakeLists.txt file discussed
# in the "CMake in user projects" page accessible from the "User info"
# page of the documentation.
SET(TARGET_SRC
  ${TARGET}.cc
  )

# Usually, you will not need to modify anything beyond this point...

CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12)

FIND_PACKAGE(deal.II 9.2.0 QUIET
  HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
  )
IF(NOT ${deal.II_FOUND})
  MESSAGE(FATAL_ERROR "\n"
"*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
"You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to 
cmake\n"
"or set an environment variable \"DEAL_II_DIR\" that contains this path."
)
ENDIF()

#
# Are all dependencies fulfilled?
#
IF(NOT ((DEAL_II_WITH_PETSC AND NOT DEAL_II_PETSC_WITH_COMPLEX) OR 
DEAL_II_WITH_TRILINOS) OR NOT DEAL_II_WITH_P4EST) # keep in one line
  MESSAGE(FATAL_ERROR "
Error! This tutorial requires a deal.II library that was configured with the 
following options:
DEAL_II_WITH_PETSC = ON
DEAL_II_PETSC_WITH_COMPLEX = OFF
DEAL_II_WITH_P4EST = ON
or
DEAL_II_WITH_TRILINOS = ON
DEAL_II_WITH_P4EST = ON
However, the deal.II library found at ${DEAL_II_PATH} was configured with these 
options
DEAL_II_WITH_PETSC = ${DEAL_II_WITH_PETSC}
DEAL_II_PETSC_WITH_COMPLEX = ${DEAL_II_PETSC_WITH_COMPLEX}
DEAL_II_WITH_P4EST = ${DEAL_II_WITH_P4EST}
DEAL_II_WITH_TRILINOS = ${DEAL_II_WITH_TRILINOS}
which conflict with the requirements.
One or both of the aforementioned combinations of prerequisites are not met by 
your installation, but at least one is required for this tutorial step."
)
ENDIF()

DEAL_II_INITIALIZE_CACHED_VARIABLES()
SET(CLEAN_UP_FILES *.log *.gmv *.gnuplot *.gpl *.eps *.pov *.vtk *.ucd *.d2 
*.vtu *.pvtu)
PROJECT(${TARGET})
DEAL_II_INVOKE_AUTOPILOT()
#include 
#include 

#include 

namespace LA {
using namespace dealii::LinearAlgebraTrilinos;
} // namespace LA

#include 
#include 
#include 

#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 

#include 
#include 
#include 
#include 

#include 
#include 
#include 

#include 
#include 

using namespace dealii;

template  class MyScalarFunction : public Function {
public:
  MyScalarFunction() : Function() {}

  virtual double value(const Point ,
   const unsigned int component = 0) const override;
  virtual void value_list(const std::vector> ,
  std::vector ,
  const unsigned int component = 0) const override;
};

template 
double MyScalarFunction::value(const Point ,
const unsigned int /*component*/) const {
  double return_value = 0;

  for (unsigned int d = 0; d < dim; ++d)
return_value += (p(d) - 0.5) * (p(d) - 0.5);

  return return_value;
}

template 
void MyScalarFunction::value_list(
const std::vector> , std::vector ,
const unsigned int /*component = 0*/) const {
  Assert(points.size() == values.size(),
 ExcDimen

[deal.II] Re: Simple merge_triangulation problem

[Konrad Wiśniewski]

Oh! Sorry for bothered you, problem lies somewhere else - this topic can be 
deleted!

-- 
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/9904bf3a-aa12-42df-8428-4fd7b6bd8bb2%40googlegroups.com.

[deal.II] Re: Simple merge_triangulation problem

[Konrad Wiśniewski]

The program is suddenly terminated with exit value -1 so, I guess it is 
segfault. It is very strange because the program works perfectly fine in 
case 1), and after I swap coordinates between cell1 and cell2 (to get the 
case 3) program just stop functioning.

W dniu czwartek, 23 stycznia 2020 15:06:48 UTC+1 użytkownik Bruno Turcksin 
napisał:
>
> Konrad,
>
> What does it mean you cannot merge the cells? Do you get the wrong 
> results? Does the code segfault? Do you get an error message?
>
> Best,
>
> Bruno
>

-- 
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/3ad806c3-ba16-4e27-81ea-e21689640501%40googlegroups.com.

[deal.II] Simple merge_triangulation problem

[Konrad Wiśniewski]

Hi all!

I've encountered a problem with very simple two-cells merging in 2D 
(deal.ii v.9.0.0). 
Lets say that I have two rectangular cells (for two different material) and 
I want to merge it together and then refine those cells properly. 
I don't have a problem with this when I want to merge cells when they are 
side by side 1):


but I cannot merge cells with this configuration 2):

Or even in this configuration 3):

I wonder why it is a case? It seems the most simple case for 
merge_triangulation() function. Can anyone show me how merge second case 
properly?

This is my code with the simplest coordinates of vertices:

//FIRST CELL
static const Point<2> vertices_1[]
= {
Point<2>(0.0, 0.0),
Point<2>(1.0, 0.0),
Point<2>(0.0, 1.0),
Point<2>(1.0, 1.0)
  };
 
// Creating CellData and vertices list for this material (in 
this case for only one cell):
const unsigned int n_vertices = 4;
const std::vector> vertices_list_2(_2[0], &
vertices_2[n_vertices]);
static const int cell_vertices[] = {0,1,2,3};
CellData cells_1;

for(unsigned int j=0;jtemp_down_triangulation;
temp_down_triangulation.create_triangulation(vertices_list_1,
cell_1,SubCellData());


//SECOND CELL
static const Point<2> vertices_2[]
= {
Point<2>(0.0, 1.0),
Point<2>(1.0, 1.0),
Point<2>(0.0, 2.0),
Point<2>(1.0, 2.0)
  };
// Creating CellData for second cell
const std::vector> vertices_list_2(_2[0], &
vertices_2[n_vertices]);
CellData cell_2;

for(unsigned int j=0;jtemp_up_triangulation;
temp_up_triangulation.create_triangulation(vertices_list_2, 
cell_2, SubCellData());

// MERGING...
GridGenerator::merge_triangulations(temp_down_triangulation,
temp_up_triangulation,
final_triangulation);


-- 
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/8286d6a4-a1f5-4191-8f81-cbfe967ebf2f%40googlegroups.com.

Re: [deal.II] Re: cmake with library and executables

[Konrad Simon]

Many thanks, Matthias!

Works!

Best,
Konrad

On Wednesday, December 25, 2019 at 2:07:41 PM UTC+1, Matthias Maier wrote:
>
>
> On Fri, Dec 20, 2019, at 13:07 CST, Konrad Simon  > wrote: 
>
> > 
> ###
>  
>
> > 
> ###
>  
>
> > ADD_CUSTOM_TARGET(debug 
> >   COMMAND ${CMAKE_COMMAND} -DCMAKE_BUILD_TYPE=Debug ${CMAKE_SOURCE_DIR} 
> >   COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target all 
> >   COMMENT "Switch CMAKE_BUILD_TYPE to Debug" 
> >   ) 
> > 
> > ADD_CUSTOM_TARGET(release 
> >   COMMAND ${CMAKE_COMMAND} -DCMAKE_BUILD_TYPE=Release 
> ${CMAKE_SOURCE_DIR} 
> >   COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target all 
> >   COMMENT "Switch CMAKE_BUILD_TYPE to Release" 
> >   ) 
> > 
> ###
>  
>
> > 
> ###
>  
>
>
> Do you get the jobserver warning while running "make release" or "make 
> debug" with "-j8"? 
>
> If so - that's because you cannot recursively call into make with the 
> jobserver feature that way. 
>
> We have fixed the example steps a while ago but probably never updated 
> the documentation regarding these two custom targets. It is best to 
> simply have: 
>
> ADD_CUSTOM_TARGET(debug 
>   COMMAND ${CMAKE_COMMAND} -DCMAKE_BUILD_TYPE=Debug ${CMAKE_SOURCE_DIR} 
>   COMMENT "Switch CMAKE_BUILD_TYPE to Debug" 
>   ) 
>
> ADD_CUSTOM_TARGET(release 
>   COMMAND ${CMAKE_COMMAND} -DCMAKE_BUILD_TYPE=Release ${CMAKE_SOURCE_DIR} 
>   ) 
>
> and then 
>
>   $ make release 
>   $ make -j8 
>
> or 
>
>   $ make debug 
>   $ make -j8 
>
> Best, 
> Matthias 
>

-- 
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/b39bba70-1033-4a18-b618-572309514c45%40googlegroups.com.

Re: [deal.II] Re: cmake with library and executables

[Konrad Simon]

Hi Wolfgang,

On Tuesday, December 24, 2019 at 5:59:31 PM UTC+1, Wolfgang Bangerth wrote:
>
>
> Konrad, 
> your email has no question :-) Is your problem that you can't call 'make 
> -j8' 
> and your question how to make that possible? If so, what happens if you 
> try? 
> What is the error message? 
>

The message essentially says I can not use the jobserver so I can only use 
once core to compile. My library does compile but falls back to make -j1 
every time I invoke make -j8. My question is why? And how can I use all 
cores?

Thank you and merry Christmas!

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/1d1a9cfa-5975-432c-9093-b542ee213cab%40googlegroups.com.

[deal.II] Re: cmake with library and executables

[Konrad Simon]

#
# Include directory for sources
#
include_directories(${MsFEComplex_INCLUDE_DIR})
###
###


###
###
#
# Name all sources
#
#file(GLOB_RECURSE MsFEComplex_TARGET_LIB_SRC  "*.cc") # source files
#set(MsFEComplex_TARGET_LIB_SRC ${MsFEComplex_TARGET_LIB_SRC})
set(MsFEComplex_TARGET_LIB_SRC
eqn_boundary_vals.cc
eqn_coeff_A.cc
eqn_coeff_B.cc
eqn_coeff_R.cc
eqn_exact_solution_lin.cc
eqn_rhs.cc
my_vector_tools.inst.cc
ned_rt_basis.cc
ned_rt_global.cc
ned_rt_post_processor.cc
ned_rt_ref.cc
parameters.cc
q_ned_basis.cc
q_ned_global.cc
q_ned_post_processor.cc
q_ned_ref.cc)

print_all_args (
  ${MsFEComplex_TARGET_LIB_SRC}
 )
###
###


###
###
#
# Compile and link the sources as SHARED
#
add_library (MsFEComplex_LIBRARY SHARED ${MsFEComplex_TARGET_LIB_SRC})
DEAL_II_SETUP_TARGET(MsFEComplex_LIBRARY)

#
# Install into the DESTINATION provided by CMAKE_INSTALL_PREFIX
#
#install (TARGETS ${MsFEComplex_LIBRARY} DESTINATION 
${CMAKE_INSTALL_PREFIX})

add_executable(MsFEComplex_Ned_RT "main_ned_rt.cxx")
DEAL_II_SETUP_TARGET(MsFEComplex_Ned_RT)
TARGET_LINK_LIBRARIES(MsFEComplex_Ned_RT MsFEComplex_LIBRARY)

add_executable(MsFEComplex_Q_Ned "main_q_ned.cxx")
DEAL_II_SETUP_TARGET(MsFEComplex_Q_Ned)
TARGET_LINK_LIBRARIES(MsFEComplex_Q_Ned MsFEComplex_LIBRARY)
###
###

So essentially I want to build a library and then link two executables to 
it. It compiles correctly but very slow since I can not invoke *make -j8*  

Google unfortunately does not help.

Any help would be much appreciated.

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/664a11a3-fcf8-46a8-92ef-42220bc8862b%40googlegroups.com.

[deal.II] projecting function onto TrilinosWrappers::MPI::BlockVector

[Konrad Simon]

Hi all,

I am having a little problem with projecting a function onto (parts of) FE 
spaces. I am getting the error

The violated condition was:  
   (dynamic_cast *>( 
&(dof.get_triangulation())) == nullptr) 
Additional information:  
   You are trying to use functionality in deal.II that is currently not 
implemented. In many cases, this indicates that there simply didn't appear 
much of a need for 
it, or that the author of the original code did not have the time to 
implement a particular case. If you hit this exception, it is therefore 
worth the time to look int
o the code to find out whether you may be able to implement the missing 
functionality. If you do, please consider providing a patch to the deal.II 
development sources 
(see the deal.II website on how to contribute).

I of course get what this error message suggests and I am wondering if I 
could fix this somehow. The funny thing is that when I step through the 
code in debug mode I see that exactly the cast above fails. Funnily, the 
cast dynamic_cast *>( 
&(dof.get_triangulation())) works. 

Now I am asking myself why?  Am I missing something here?

Best regards,
Konrad


This is my function:

{
TrilinosWrappers::MPI::BlockVector locally_relevant_exact_solution;
locally_relevant_exact_solution.reinit(owned_partitioning,
mpi_communicator);

{ // write sigma to Nedelec_0 space

// Quadrature used for projection
QGauss<3> quad_rule (3);

// Setup function
ExactSolutionLin_A_curl exact_sigma(parameter_filename); // Exact solution 
for first FE ---> Nedelec FE

DoFHandler<3> dof_handler_fake (triangulation);
dof_handler_fake.distribute_dofs (fe.base_element(0));

if (parameters.renumber_dofs)
{
DoFRenumbering::Cuthill_McKee (dof_handler_fake);
}

AffineConstraints constraints_fake;
constraints_fake.clear ();
DoFTools::make_hanging_node_constraints (dof_handler_fake, 
constraints_fake);
constraints_fake.close();

VectorTools::project (dof_handler_fake,
constraints_fake,
quad_rule,
exact_sigma,
locally_relevant_exact_solution.block(0));

dof_handler_fake.clear ();
}

{// write u to Raviart-Thomas_0 space

// Quadrature used for projection
QGauss<3> quad_rule (3);

// Setup function
ExactSolutionLin exact_u(parameter_filename); // Exact solution for second 
FE ---> Raviart-Thomas FE

DoFHandler<3> dof_handler_fake (triangulation);
dof_handler_fake.distribute_dofs (fe.base_element(1)); 

if (parameters.renumber_dofs)
{
DoFRenumbering::Cuthill_McKee (dof_handler_fake);
}

AffineConstraints constraints_fake;
constraints_fake.clear ();
DoFTools::make_hanging_node_constraints (dof_handler_fake, 
constraints_fake);
constraints_fake.close();

VectorTools::project (dof_handler_fake,
constraints_fake,
quad_rule,
exact_u,
locally_relevant_exact_solution.block(1));

dof_handler_fake.clear ();
}
}

-- 
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/7a66b73f-4fe3-449f-b0ba-c4b16e8fbca3%40googlegroups.com.

[deal.II] Question about postprocessing

[Konrad Simon]

Hi deal.ii community,

I also have a little question about postprocessing in different spaces. I 
am post processing two solutions of the same problem but solved in two 
different (pairs of) spaces. One quantity, for example, is called u and is 
either in H(curl) or in H(div) depending on the form of the problem. 
Another is sigma and it is either in H^1 or in H(curl), respectively.

I need to compute things that involve the divergence of u when u is in 
H(div) or the curl when u is in H(curl). Both things I do using the 
DataPostprocessor class with entries of the (matrix valued) gradient of u 
but I don't know the internals.

My problem ist that when comparing quantities that should be similar 
according to the math then I get differences that are too large 
(intuitively). 

The thing is that when I have a quantity that is in H(curl) using Nedelec 
approximation then I can nicely take the curl but divergence will be zero 
by construction. Same if u is in H(div) with Raviart-Thomas approximation 
(then the curl is zero by construction). How is the computation of the 
gradient done internally?

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/68c66960-82d8-4a63-92ac-9c104926959a%40googlegroups.com.

[deal.II] cmake with library and executables

[Konrad Simon]

Hi deal.ii community,

Little cmake question: I set up a user project with include, test, doc and 
source directory. I am collecting some code in a library and then I link a 
few application files to the library. I followed the guidelines of the 
documentation and everything is fine. Now one little but annoying thing. I 
get

make[4]: warning: jobserver unavailable: using -j1.  Add '+' to parent make 
rule.

when typing make -j8.

Do you have an idea what I could do?

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/8d5a7e5e-3cc1-4fb9-bd3d-63f6db4622f8%40googlegroups.com.

Re: [deal.II] Evaluate shape functions for an element on a given cell in a triangulation

[Konrad Simon]

Many thanks, Wolfgang. 

-- 
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/8ce85110-2bf7-4f76-a4cb-ac1f9bdf3781%40googlegroups.com.

[deal.II] Re: How to use CellId class

[Konrad Simon]

Hi Zhidong,

On Monday, October 21, 2019 at 1:42:30 AM UTC+2, Zhidong Brian Zhang wrote:
>
> Thank you very much for your prompt reply, Konrad!
>
> My confusion is the output of cell->id(), for example,
> 0_3:000
> 0_3:200
> 0_3:003
> 0_3:006
> 0_3:406
> 0_3:606
> 0_3:206
> 0_3:007
> 0_3:407
> 0_3:607
> 0_3:207.
>
> What type I need to use as the key in std::map? And is the map 
> parallel-distributed? Because my model requires that.
>
> You can use the CellId object as a key itself: 

std::map cell_info_map;

If you need the map distributed across nodes according to the distribution 
of cells then you can initialize the map (on each node separately) through

for (; cell!=endc; ++cell)
{
if (cell->is_locally_owned())
{
CellId current_cell_id(cell->id());

std::pair::iterator, bool > 
result;
result = cell_basis_map.insert(std::make_pair(cell->id(),
cell_info_map));

 // ClassWithCellInfo must have a copy constructor
Assert(result.second,
ExcMessage ("Insertion of cell_info_map into std::map failed. "
"Problem with copy constructor?"));
 } // end if
} // end ++cell

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/51256993-cd1b-487d-82a4-ee08ce56c501%40googlegroups.com.

[deal.II] Evaluate shape functions for an element on a given cell in a triangulation

[Konrad Simon]

Hi deal.ii community,

I have a little implementation problem to bug you with... 

Is there a way to evaluate a given shape function on a given cell like a 
normal Function or TensorFunction? I need to do many computations on a 
coarse cell that itself is meshed.

I have an implementation but it is not very efficient: I create an FEValues 
object and first map the evaluation points back to the unit cell using an 
appropriate mapping. These I use as "quadrature points" and get the values 
on the real (coarse) cell by then reinitializing the FEValues object with 
it. Works for all shape functions but is slow and cumbersome.

I actually only need it for Q1 elements and lowest order Nedelec and 
Raviart-Thomas elements.

Anyone has an idea?

Cheers,
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/c3226e74-c327-4f63-8be9-ba71468a3a7e%40googlegroups.com.

[deal.II] Re: How to use CellId class

[Konrad Simon]

Hi Zhidong,

I don't know what exactly you want to do but I found it quite useful that 
the CellId class has an order relation. That makes it usable in a std::map. 
Use the CellId (you can retrieve it through cell-->id() if cell is a cell 
accessor like the one you use to loop over all your cells) as a key then 
and connect to it an object that, for example, contains specific 
information about your cell(s).

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/829003a9-f318-4790-b957-9684aec7b91e%40googlegroups.com.

[deal.II] Re: Issue with boost serialization and spack?

[Konrad Simon]

Hi Denis,

I don't have the build folders any more so I can not post the error log. 
But the error (using spack) occurred with both boost versions. I will post 
something once I will find a solution.
At any rate, thank you for your help.

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/4fab9bae-8ed2-49de-918c-a270a4070c33%40googlegroups.com.

[deal.II] Re: Issue with boost serialization and spack?

[Konrad Simon]

>
>
>> I guess I have a clue why I get the error. My backend nodes run on a 
>> different architecture. The openmpi version (i.e. mpirun) is the correct 
>> one. This is made sure in the slurm script.
>>
> That's always a pain to deal with. If you can, I would get an interactive 
> session on a compute node and compile everything on the compute node. You 
> can also compile everything in batch mode or compute on the login node and 
> pass the flag for the architecture of the compute node.
>
> Thank you, that's what I am trying. I will find a solution the question is 
jsut how much time it will eat...
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/f3c230bb-68b1-4491-ae3b-04ec9f6c74f8%40googlegroups.com.

[deal.II] Re: Issue with boost serialization and spack?

[Konrad Simon]

>
> Thank you, Denis. I use a pretty stupid (but simple) workaround: I setup 
>> and compile deal.ii myself since all dependencies are compiled and use the 
>> cmake command used by spack. That works. And I do not get the serialization 
>> error.
>>
>
>
> now that is strange. Are you sure you pick up the same boost?
> Could you post the original error from CMake error logs that shows how 
> serialization test fails?
>  
>
Cmake in my little modification seems to pick up the right boost. (for 
error log with boost see above)
 
DEAL_II_WITH_BOOST set up with external dependencies 
#BOOST_VERSION = 1.70.0 
#BOOST_DIR = 
/scratch/cen/numgeo/spack-lib/linux-debian8-x86_64/gcc-7.3.0/boost-1.70.0 
#BOOST_CXX_FLAGS = -Wno-unused-local-typedefs 
#BOOST_DEFINITIONS = BOOST_NO_AUTO_PTR 
#BOOST_USER_DEFINITIONS = BOOST_NO_AUTO_PTR 
#BOOST_INCLUDE_DIRS = 
/scratch/cen/numgeo/spack-lib/linux-debian8-x86_64/gcc-7.3.0/boost-1.70.0/include
 

#BOOST_USER_INCLUDE_DIRS = 
/scratch/cen/numgeo/spack-lib/linux-debian8-x86_64/gcc-7.3.0/boost-1.70.0/include
 

#BOOST_LIBRARIES = 
/scratch/cen/numgeo/spack-lib/linux-debian8-x86_64/gcc-7.3.0/boost-1.70.0/lib/libboost_iostreams-mt.so;/scratch/cen/numgeo/spack-lib/lin
ux-debian8-x86_64/gcc-7.3.0/boost-1.


 

>
>> However, now my code runs on the machine I installed it on. But once I 
>> use slurm to distribute the job across nodes I get "illegal instruction" 
>> erros. Frustrating.
>>
>
> With HPC I used to have access to they did not use slurm, so I can't 
> really comment here.
> If that's during running quick tests or so, it could be related to a wrong 
> MPIEXEC pickedup by deal.II config.
> I wanted to fix that in Spack https://github.com/spack/spack/pull/11142 
> but apparently this solution may not be fully functional for Slurm.
>

I guess I have a clue why I get the error. My backend nodes run on a 
different architecture. The openmpi version (i.e. mpirun) is the correct 
one. This is made sure in the slurm script.
 
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/d1e8e796-969e-4c41-ae0c-1c5a54742dd2%40googlegroups.com.

[deal.II] Data postprocessing with two dof handlers

[Konrad Wiśniewski]

 

 Hi! 


I'm looking for quick advise. In my calculations I solve two equations (in 
2D) consecutively: the Poisson equation where I obtain electric field *E*, 
and then I solve the continuity equation where I obtain e.g. density of 
electrons *n*. I want to get current that is *j* = C*nE* where C is a 
constant (I need to use one variable from first and one from second 
equation). Problem is that in Poisson equation I've used Raviart-Thomas 
finite elements and in Continuity equation I've used discontinuous elements 
(FE_DGQ) and the position of dofs in this two methods do not coincide with 
each other (RT dofs are placed on the faces of a cell when DGQ on the 
vertexes). My question is: how do you wisely combine two dof handlers to 
produce a new variable which I want to visualize on the same mesh?

E.g: How to extrapolate the values from RT-elements to vertexes, then for 
each vertex calculate current *j*, and then add current variable to the 
solution in the way that I will be able to visualize it in paraview?


Bonus question: I'm not sure now but I may need to have two different 
triangulations in those equations, so even the cells will differ one from 
another (however the overall domain will remain the same). How should I 
deal with this situation?

-- 
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/595dcf90-2713-427a-ab39-d055c4ed6f7d%40googlegroups.com.

[deal.II] Re: Issue with boost serialization and spack?

[Konrad Simon]

Thank you, Denis. I use a pretty stupid (but simple) workaround: I setup 
and compile deal.ii myself since all dependencies are compiled and use the 
cmake command used by spack. That works. And I do not get the serialization 
error.

However, now my code runs on the machine I installed it on. But once I use 
slurm to distribute the job across nodes I get "illegal instruction" erros. 
Frustrating.

Nevertheless, thanks for your hint.

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/49a66335-e58d-4eab-aee3-55e87a54f88f%40googlegroups.com.

[deal.II] Issue with boost serialization and spack?

[Konrad Simon]

Dear deal.ii community,

I am having a little problem and I was wondering if this issue is known. I 
installed deal.ii on our cluster  and all dependencies build nicely with my 
chosen compiler (gcc v7.3). BLAS and LAPACK are being built as well as 
openmpi (the versions on the cluster are not suitable at the moment). 
However, when building deal.ii I get the output below (I only paste the 
relevant part). It also happens when I ask preck to build deal.ii with 
(externally spack built) boost v1.62 and v1.68. 

Could there be a configuration issue with boost in the spack build class of 
deal.ii? When looking at it I see that some version dependent boost patches 
are applied but I am not familiar with the details.

Anyone knows this? Is there a workaround?

Thanks and best regards,
Konrad

-- Include 
/tmp/u290231/spack-stage/dealii-9.1.1-gfk33pgt5rojhujhmuruaahuuzyzq2zm/spack-src/cmake/configure/configure_2_boost.cmake
-- Found Boost: 
/scratch/cen/numgeo/spack-lib/linux-debian8-broadwell/gcc-7.3.0/boost-1.62.0/include
 
(found suitable version "1.62.0", minimum required is "1.59") found 
components:  iostreams serialization system thread regex chrono date_time 
atomic 
--   BOOST_VERSION: 1.62.0
--   BOOST_LIBRARIES: 
/scratch/cen/numgeo/spack-lib/linux-debian8-broadwell/gcc-7.3.0/boost-1.62.0/lib/libboost_iostreams-mt.so;/scratch/cen/numgeo/spack-lib/linux-debian8-broadwell/gcc-7.3.0/boost-1.62.0/lib/libboost_serialization-mt.so;/scratch/cen/numgeo/spack-lib/linux-debian8-broadwell/gcc-7.3.0/boost-1.62.0/lib/libboost_system-mt.so;/scratch/cen/numgeo/spack-lib/linux-debian8-broadwell/gcc-7.3.0/boost-1.62.0/lib/libboost_thread-mt.so;/scratch/cen/numgeo/spack-lib/linux-debian8-broadwell/gcc-7.3.0/boost-1.62.0/lib/libboost_regex-mt.so;/scratch/cen/numgeo/spack-lib/linux-debian8-broadwell/gcc-7.3.0/boost-1.62.0/lib/libboost_chrono-mt.so;/scratch/cen/numgeo/spack-lib/linux-debian8-broadwell/gcc-7.3.0/boost-1.62.0/lib/libboost_date_time-mt.so;/scratch/cen/numgeo/spack-lib/linux-debian8-broadwell/gcc-7.3.0/boost-1.62.0/lib/libboost_atomic-mt.so
--   BOOST_INCLUDE_DIRS: 
/scratch/cen/numgeo/spack-lib/linux-debian8-broadwell/gcc-7.3.0/boost-1.62.0/include
--   BOOST_USER_INCLUDE_DIRS: 
/scratch/cen/numgeo/spack-lib/linux-debian8-broadwell/gcc-7.3.0/boost-1.62.0/include
-- Found BOOST
-- Performing Test BOOST_IOSTREAMS_USABLE
-- Performing Test BOOST_IOSTREAMS_USABLE - Success
-- Performing Test BOOST_SERIALIZATION_USABLE
-- Performing Test BOOST_SERIALIZATION_USABLE - Failed
-- The externally provided Boost.Serialization library failed to pass a 
crucial test. 
Therefore, the bundled boost package is used. 
The configured testing project can be found at 
/tmp/u290231/spack-stage/dealii-9.1.1-gfk33pgt5rojhujhmuruaahuuzyzq2zm/spack-build/cmake/configure/TestBoostBugWorkdir
-- DEAL_II_WITH_BOOST has unmet external dependencies.
CMake Error at cmake/macros/macro_configure_feature.cmake:112 (MESSAGE):
  

  Could not find the boost library!

  The externally provided Boost.Serialization library failed to pass a
  crucial test.Please ensure that a suitable boost library is installed on
  your computer.

  If the library is not at a default location, either provide some hints for
  autodetection,

  $ BOOST_DIR="..." cmake <...>
  $ cmake -DBOOST_DIR="..." <...>

  or set the relevant variables by hand in ccmake.

  Alternatively you may choose to compile the bundled library of boost by
  setting DEAL_II_ALLOW_BUNDLED=on or DEAL_II_FORCE_BUNDLED_BOOST=on.

Call Stack (most recent call first):
  cmake/macros/macro_configure_feature.cmake:269 (FEATURE_ERROR_MESSAGE)
  cmake/configure/configure_2_boost.cmake:237 (CONFIGURE_FEATURE)
  cmake/macros/macro_verbose_include.cmake:19 (INCLUDE)
  CMakeLists.txt:124 (VERBOSE_INCLUDE)

-- 
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/e889b864-66ad-4e05-81f3-08420d5e6a30%40googlegroups.com.

[deal.II] Re: Configuration of PETSc

[Konrad Simon]

Hi Toni,

Seems like I missed that little note in the documentation. Thank you :-)

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/3148cde3-83f1-46a4-be88-3ddfe3e45cdc%40googlegroups.com.

[deal.II] Re: deal.II for Inverse problems in non-linear elasticity

[Konrad Simon]

Hi Prashant,
 

> I am trying to solve an Inverse Cauchy problem in 3D nonlinear elasticity. 
> I have observed displacement data at partial boundary as well in partial 
> regions inside the body, and want to reconstruct the traction field. From 
> documentation of deal.II, I understood how tangent stiffness matrix can be 
> created at each iteration for the forward problem.
>
 
This is a challenging ill-posed problem. I did something similar myself in 
the past. An advice from my experience: Start with linear elasticity 
(Navier-Lame). In fact the stiffness-matrix is the tangent stiffness matrix 
of Saint-Venant Materials around the zero-displacement field. Then with 
Saint-Venant. If this works try "simple" Neo-Hookean laws (keep in mind 
that tangent stiffness is not necessarily definite etc). etc... but I guess 
that was not your question. (Sorry for the unrequested advice)

While solving the inverse problem iteratively,  I want to retrieve this 
> tangent stiffness matrix at each iteration, perform necessary algebraic 
> manipulations and update the increments in unknown displacement and 
> traction variables.  Is it possible to do so using deal.II?  If you have 
> any suggestions, please feel free to include them.
>

Look at tutorial step-44. I think you will find what you are looking for.

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/56f9a37a-58e0-4b56-a265-41ebdd1557b0%40googlegroups.com.

[deal.II] Configuration of PETSc

[Konrad Simon]

Dear deal.ii community,

I am using deali.ii with PETSc and Trilinos. However, when I am using the 
PETSc PreconditionILU I get an error that suggests that a solver package is 
missing (with Trilinos it works). Petsc's PreconditionAMG works fine 
(although not very efficiently for my problem). 
Do I need to do any special configuration steps for PETSC? I followed the 
instructions that are documented on the deal.ii pages on "how to configure 
Petsc". https://www.dealii.org/current/external-libs/petsc.html

Best,
Konrad

This is the error:

Running using PETSc. 
Number of active cells: 262144 
Total number of cells: 24865 (on 7 levels) 
Number of degrees of freedom: 1609920 (811200+798720) 
[0]PETSC ERROR: - Error Message 
-- 
[0]PETSC ERROR: See 
http://www.mcs.anl.gov/petsc/documentation/linearsolvertable.html for 
possible LU and Cholesky solvers 
[0]PETSC ERROR: Could not locate a solver package. Perhaps you must 
./configure with --download- 
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for 
trouble shooting. 
[0]PETSC ERROR: Petsc Release Version 3.9.4, Sep, 11, 2018  
[0]PETSC ERROR: main on a x86_64 named thunder5 by u290231 Fri Sep 20 
10:00:23 2019 
[0]PETSC ERROR: Configure options --with-shared-libraries=1 --with-x=0 
--with-mpi=yes --download-hypre=yes --with-64-bit-indices 
--with-debugging=yes --with-hypre=yes 
[0]PETSC ERROR: #1 MatGetFactor() line 4328 in 
/scratch/cen/numgeo/lib/petsc-3.9.4/src/mat/interface/matrix.c 
[0]PETSC ERROR: #2 PCSetUp_ILU() line 142 in 
/scratch/cen/numgeo/lib/petsc-3.9.4/src/ksp/pc/impls/factor/ilu/ilu.c 
[0]PETSC ERROR: #3 PCSetUp() line 923 in 
/scratch/cen/numgeo/lib/petsc-3.9.4/src/ksp/pc/interface/precon.c 
- 
TimerOutput objects finalize timed values printed to the 
screen by communicating over MPI in their destructors. 
Since an exception is currently uncaught, this 
synchronization (and subsequent output) will be skipped to 
avoid a possible deadlock. 
- 
WARNING! There are options you set that were not used! 
WARNING! could be spelling mistake, etc! 
Option left: name:-p value: ../MsFEComplex/parameter.in 
ERROR: Uncaught exception in MPI_InitFinalize on proc 0. Skipping 
MPI_Finalize() to avoid a deadlock. 


 
Exception on processing:  

 
An error occurred in line <421> of file 
 
in function 
   void dealii::PETScWrappers::PreconditionILU::initialize(const 
dealii::PETScWrappers::MatrixBase&, const 
dealii::PETScWrappers::PreconditionILU::AdditionalData&) 
The violated condition was:  
   ierr == 0 
Additional information:  
deal.II encountered an error while calling a PETSc function. 
The description of the error provided by PETSc is "See 
http://www.mcs.anl.gov/petsc/documentation/linearsolvertable.html for 
possible LU and Cholesky solvers". 
The numerical value of the original error code is 92. 
 

Aborting! 
 
-- 
Primary job  terminated normally, but 1 process returned 
a non-zero exit code. Per user-direction, the job has been aborted. 
-- 
-- 
mpirun detected that one or more processes exited with non-zero status, 
thus causing 
the job to be terminated. The first process to do so was: 

 Process name: [[37149,1],0] 
 Exit code:1 
--

-- 
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/71e73e1e-0a03-417d-b086-89b62e44a2e1%40googlegroups.com.

Re: [deal.II] Initialization of iterative solver

[Konrad Simon]

>
> Yes, you just set the solution vector to your best guess before you pass 
> it to the solver. But you'll find that in practice, this rarely reduced 
> the number of iterations by any significant amount :-( 
>

Many thanks, Wolfgang. Now, I also saw that in the source code of 
solver_cg.h. I really hoped that would speed up things a bit when solving 
time dependent problems.. Hm. 

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/c1d3ff11-38c3-4c03-b8a0-5329b3247c27%40googlegroups.com.

[deal.II] Initialization of iterative solver

[Konrad Simon]

Dear deal.ii community,

I have a quick question: Suppose I know a vector that is close to the 
solution of my linear system. Is there a method to initialize an iterative 
solver with this vector?

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/c9a7e762-bf1b-4ce3-91c1-515610ffdb5b%40googlegroups.com.

[deal.II] Re: On using project_boundary_values_div_conforming for (RT-DG) in parallel case

[Konrad]

Hi Charlie, 

I can compile and run it. Could there be a problem with your deal.ii setup? 
I am using v9.1.1. 

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/47f7594a-e5e5-44dd-99c0-1c41a1757992%40googlegroups.com.

Re: [deal.II] Question about (parallel) linear solver

[Konrad]

Hi Wolfgang,
 

> > Using approach 1.) I get timing results like that: 
> > (using MPI with Petsc, AMG preconditioning for *A* on one machine with 
> > four cores, no preconditioning for the Schur complement *S = 
> > -B*A^{-1}*B^T + C* from 1.) 
> > 
> > Running using PETSc. 
> > Number of active cells: 32768 
> > Total number of cells: 21449 (on 6 levels) 
> > Number of degrees of freedom: 205920 (104544+101376) 
> >Iterative Schur complement solver converged in 177 iterations. 
> >Outer solver completed. 
> > 
> > 
> > 
> +-+++ 
> > | Total wallclock time elapsed since start|37s | 
>| 
> > | || 
>| 
> > | Section | no. calls |  wall time | % of total 
> | 
> > 
> +-+---+++ 
> > | Schur complement solver (for u) | 1 |  27.2s |74% 
> | 
> > | assembly| 1 |  6.61s |18% 
> | 
>
> That's actually an acceptable time: My rule of thumb, maybe a bit 
> outdated, is that it takes a minute to solve for 100,000 unknowns for 
> non-trivial problems. You've got twice that many unknowns and solve it 
> in half the time. 
>
Thank you, that calms me down.
 

> > You will notice that the time for the Schur complement solver is 
> > dominating. So either I did something wrong in the implementation 
> > conceptually (see class for Schur complement below) or I need to do a 
> > better job in solving/preconditioning the system. Note that the Schur 
> > complement is usually not definite in contrast to Stokes or Darcy 
> systems. 
>
> Really? Is this an issue related to boundary conditions? Your matrix 
> corresponds to the operator 
>curl I curl + grad div 
> I thought there are no vector fields that are both divergence-free and 
> curl-free. 
>

That indeed depends on the geometry of the domain and it is not an obvious 
result. Example: Think of a 3D hyper-shell with an inner and outer radius. 
There is a one-dimensional space of non-trivial radial vector fields (with 
respect to the origin) that do not have neither a divergence nor a curl. 
Such fields are referred to as harmonic forms and the dimensionality of 
this space depends on the Betti numbers of the domain (buzz word is Hodge 
decomposition).

But this is actually not the case in here, so you are right here. If there 
are non-trivial harmonic forms one must have an additional condition and a 
Lagrange multiplier in the weak form.
 

> > Does anyone have an idea of how to attack systems like *(P)*? Are there 
> > optimal preconditioners? Any hint would be appreciated. 
>
> Take a look at step-20 and step-22. The idea of preconditioning a Schur 
> complement is to use some sort of approximation of the Schur complement 
> that is easier to solve for. I suspect that in your case, a simple 
> Laplace matrix would make for a nice approximation of S, and so using 
> one SSOR step with the Laplace matrix might be a good preconditioner for 
> S. Alternatively, an algebraic multigrid initialized with the Laplace 
> matrix might work. 
>
> step-20 plays with these ideas, solving for the Schur complement 
> directly. step-22 (in particular the code in the "Possibilities for 
> extensions" section) expands on this: If you have a preconditioner for 
> the Schur complement, then you can just use that as part of the 
> preconditioner for the overall system. It's the best approach we have 
> for the Stokes equation (used successfully on systems with billions of 
> unknowns) and worth understanding! 
>

This is great to know and I will then go in this direction. Again, many 
thanks :-)

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/e92249cc-9464-436e-97a2-36bf3710ebca%40googlegroups.com.

[deal.II] Re: On the usage of Utilities::MPI

[Konrad]

>
> However, when I run the command "mpirun -np 5 ./MPI_TEST", I get the 
> following message,
>
> Hello world from process 0 of 1
> Hello world from process 0 of 1
> Hello world from process 0 of 1
> Hello world from process 0 of 1
> Hello world from process 0 of 1
>
> It seems that 5 cpu cores are used, however, they seemed to have the same 
> rank. 
> Besides, my operating system is Ubuntu 16.04. Looking forward to hearing 
> from you soon. Thank you!
>

Hi Yang Liu,

Try that code (explanations below): 

#include 
#include 
#include 

using namespace dealii;

template 
class MPI_TEST
{
public:
  MPI_TEST();
  void print_test();

private:
  MPI_Comm mpi_communicator;
  ConditionalOStream pcout;  
};

template 
void MPI_TEST::print_test()
{
  for (unsigned int i = 0;
  i < Utilities::MPI::n_mpi_processes(mpi_communicator);
  ++i)
  {
  std::cout << "Hello world from process "
  << Utilities::MPI::this_mpi_process(mpi_communicator)
  << " of "
  <<  Utilities::MPI::n_mpi_processes(mpi_communicator)
  << std::endl;
  }
}

template 
MPI_TEST::MPI_TEST()
  : mpi_communicator(MPI_COMM_WORLD)
  , pcout(std::cout, (Utilities::MPI::this_mpi_process(mpi_communicator) == 
0))
{}


int main(int argc, char ** argv)
{

  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);

  MPI_TEST<2> hw_mpi;
  hw_mpi.print_test();

  return 0;
}

The out put of mpirun -np 3 ./MPI_TEST is:

Hello world from process 1 of 3
Hello world from process 1 of 3
Hello world from process 1 of 3
Hello world from process 0 of 3
Hello world from process 0 of 3
Hello world from process 0 of 3
Hello world from process 2 of 3
Hello world from process 2 of 3
Hello world from process 2 of 3

Looks weird if you are not used to MPI but it is as expected. 

1. You must call Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, 
argv, 1); before you invoke any action of MPI. It will create an object 
that calls MPI_Init upon creation and MPI_Finalize in the destructor, i.e., 
when the object goes out of scope. Only then

2. When you implement MPI code it is important to get into the right 
mindset: You do not program code for one compute node only. You are 
programming code for all nodes (at the same time). In your example the for 
loop gets executed on all nodes( =3 here).

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/5164476e-19e6-4e08-a0a2-f938a9adab89%40googlegroups.com.

[deal.II] Question about (parallel) linear solver

[Konrad]

Dear deal.ii community,

Apologies for so many questions lately... but here is another one. This one 
is a bit technical and aims at good solvers. In a simplified form you can 
formulate it this way:

I implemented a PDE solver with a for the 3D vector Laplacian *curl(curl u) 
- grad(div u) = f*. I use a Nedelec-Raviart-Thomas pairing to solve the 
mixed form which reads
*sigma - curl u = 0*
*curl sigma - grad(div u) = f*

You can think of this as a Helmholtz decomposition of the vector field *f*. 
I am using Nedelec elements for *sigma* and Raviart-Thomas elements for *u*. 
Fine so far. Besides the fact that this is technically not so easy when it 
comes to certain domains with holes etc (forget about non-trivial harmonic 
forms if this tells you something - I don't have that here) it seems quite 
challenging to come up with a good solver (what is "good" may depend on *f* 
as well).

The system that needs to be solved formally reads

*(P)   A sigma + B^T u = 0*
*B sigma + C u = f*

*A* is sort of a mass matrix so it is positive (corresponds to ** 
in *H(curl)*) and *C* is positive semi-definite (corresponds to ** in *H(div)*).

First thing that crossed my mind was a Schur complement solver. Now one a 
priori has two options:

1.) Solve first for u and then for sigma, i.e. invert *A* and use the Schur 
complement *S = -B*A^{-1}*B^T + C*
2.) The other way around *S = -B^T*C^{-1}*B + A*. My gut feeling tells me 
that this could actually turn out to be not an option since either *C* is 
very ill-conditioned or not invertible at all. It is symmetric and contains 
at least one zero eigenvalue. Thinking of it: it formally corresponds to 
*C~(grad 
div)^{-1}* (and I am of course having convergence issues)

Using approach 1.) I get timing results like that:
(using MPI with Petsc, AMG preconditioning for *A* on one machine with four 
cores, no preconditioning for the Schur complement *S = -B*A^{-1}*B^T + C* 
from 1.)

Running using PETSc. 
Number of active cells: 32768 
Total number of cells: 21449 (on 6 levels) 
Number of degrees of freedom: 205920 (104544+101376) 
  Iterative Schur complement solver converged in 177 iterations. 
  Outer solver completed. 


+-+++ 
| Total wallclock time elapsed since start|37s || 
| ||| 
| Section | no. calls |  wall time | % of total | 
+-+---+++ 
| Schur complement solver (for u) | 1 |  27.2s |74% | 
| assembly| 1 |  6.61s |18% | 
| mesh generation | 1 | 0.163s |  0.44% | 
| outer CG solver (for sigma) | 1 | 0.136s |  0.37% | 
| system and constraint setup | 1 |  1.15s |   3.1% | 
| vtu output  | 1 |  1.68s |   4.5% | 
+-+---+++ 


You will notice that the time for the Schur complement solver is 
dominating. So either I did something wrong in the implementation 
conceptually (see class for Schur complement below) or I need to do a 
better job in solving/preconditioning the system. Note that the Schur 
complement is usually not definite in contrast to Stokes or Darcy systems.

Does anyone have an idea of how to attack systems like *(P)*? Are there 
optimal preconditioners? Any hint would be appreciated. Are there relevant 
tutorials? Parallelization is an issue since I would like to solve very 
large problems in 3D.

Thanks in advance and best regards,
Konrad :-)

Here the Schur complement code:

template 
class SchurComplementMPI : public Subscriptor
{
private:
using BlockType = typename BlockMatrixType::BlockType;

public:
SchurComplementMPI (const BlockMatrixType _matrix,
   const InverseMatrix 
_inverse_matrix,
   const std::vector _partitioning,
   MPI_Comm mpi_communicator);

void vmult (VectorType   ,
  const VectorType ) const;

private:
const SmartPointer system_matrix;
const SmartPointer> 
relevant_inverse_matrix;

const std::vector& owned_partitioning;

MPI_Comm mpi_communicator;

mutable VectorType tmp1, tmp2, tmp3;
};


template 
SchurComplementMPI::SchurComplementMPI(
const BlockMatrixType _matrix,
const InverseMatrix _inverse_matrix,
const std::vector _partitioning,
MPI_Comm mpi_communicator)
:
system_matrix (_matrix),
relevant_inverse_matrix (_inverse_matrix),
owned_partitioning(owned_partitioning),
mpi_communicator(mpi_communicator),
tmp1 (owned_partitioning[0],
mpi_communicator),
tmp2 (owned_partitioning[0],
mpi_communicator),
tmp3 (owned_partitioning[1],
mpi_communicator)
{}


template 
void
SchurComplementMPI::vmult(
VectorType   ,
const VectorType ) const
{
system_matrix->block(0,1).vmult (tmp1, src);
relevant_inverse_matrix->vmult (tmp2, tmp1

Re: [deal.II] Re: Error with Petsc and and multithreading?

[Konrad]

>
> Yes, this is indeed the case, and is what the error message was trying to 
> suggest. If you think that the error message could have been clearer, feel 
> free to propose a better wording and we'll include it for the next 
> release! 
>

Thanks a lot, Wolfgang! The error message is actually quite clear. What I 
found out in addition is that it is advantageous and helpful to read the 
documentation ... ;-)

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/6aa4d4eb-f91c-4519-a707-4833225fe4d8%40googlegroups.com.

[deal.II] Re: Error with Petsc and and multithreading?

[Konrad]

OK, just to answer the question: If you are running with petsc you should 
not call

dealii::Utilities::MPI::MPI_InitFinalize
mpi_initialization(argc, argv, dealii::numbers::invalid_unsigned_int);

since this will invoke threading if you start a number of mpi processes 
that is not equal the number of cores you use. Instead one should pass 
(when using generic linear algebra)

#ifdef USE_PETSC_LA
dealii::Utilities::MPI::MPI_InitFinalize
mpi_initialization(argc, argv, /* disable threading for petsc */ 1);
#else
dealii::Utilities::MPI::MPI_InitFinalize
mpi_initialization(argc, argv, dealii::numbers::invalid_unsigned_int);
#endif

I guess that was the bug.

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/628d94ee-9e38-4dfd-b0b7-0a84d7b832ac%40googlegroups.com.

Re: [deal.II] Poisson equation, implementation of boundary values, mixed finite element, semiconductor devices

[Konrad Wiśniewski]

Hi Konrad,

Thank you very much for your comment! I will try to go into the detail of 
it in next days. I must admit that to this point I've just used 
::make_zero_boundary_constraints with regard to the conditions on E field. 
I didn't have a better idea, and your solutions seem much more better. I 
will probably use the second option, but I need to read more about the 
first function that you mentioned. Anyway - it will help me for sure, 
thanks! We will see if I will be able to proceed without any further advice 
:). 

Best
Konrad



W dniu sobota, 7 września 2019 09:42:10 UTC+2 użytkownik Konrad napisał:
>
> Hi Konrad,
>
> I may add as a comment (without going into details): Your problem is - as 
> you correctly describe above - a mixed Poisson problem. You essentially 
> have two options to solve it and the way you choose influences the way you 
> need to treat the boundary conditions: in primal form and in mixed form. In 
> your formulation you have natural boundary conditions (Dirichlet) on Phi 
> which means that you need to prescribe them in the variational form. This 
> will lead to a boundary integral on the relevant part of the boundary and 
> result in an additional term on the right-hand side. The boundary condition 
> for E is a flux condition and in your mixed form this (Neumann) condition 
> is essential. This means that you need to enforce them on the 
> Raviart-Thomas space as a hard constraints. There are essentially two 
> options for this:
>
> 1. If you have given a vector field that describes E on the relevant part 
> of the boundary you can use the function
> *VectorTools::project_boundary_values_div_conforming. *This will project 
> the given vector field onto its normal component (since this is the one you 
> must describe with RT elements in H(div)) and use the projected value as 
> boundary condition for the normal flux of E
>
> 2. If you only are given the scalar value of the normal flux of E then you 
> may have to set the relevant dofs in the RT-space by hand. This you can do 
> through constraints that way:
>
> AffineConstraints   constraints;
>
> constraints.clear ();
>
> DoFTools::make_hanging_node_constraints (dof_handler, constraints);
>
> const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
> const unsigned int   dofs_per_face   = fe.dofs_per_face;
> std::vector local_dof_indices (dofs_per_cell);
> std::vector local_dof_face_indices 
> (dofs_per_face);
>
> typename DoFHandler::active_cell_iterator
> cell = dof_handler.begin_active(),
> endc = dof_handler.end();
> for (; cell!=endc; ++cell)
> {
> cell->get_dof_indices(local_dof_indices);
>
> for (unsigned int face_n=0;
> face_n::faces_per_cell;
> ++face_n)
> {
> cell->face(face_n)->get_dof_indices(local_dof_face_indices);
> if (cell->at_boundary(face_n))
> {
> if (cell->face(face_n)->boundary_id()==id_bottom_boundary)
> {
> for (unsigned int i = 0; i {
> // see documentation for these functions
> constraints.add_line (local_dof_face_indices.at
> (i));
> constraints.set_inhomogeneity (
> local_dof_face_indices.at(i),
> value_of_dof);
> }
> }
> else // add some constraints on other boundaries. Here: 
> zero-flux
> {
> for (unsigned int i = 0; i {
> constraints.add_line (local_dof_face_indices.at
> (i));
> }
> }
> }
> }
> ++cell_n;
> }
>
> constraints.close();
>
>
> Once you use the DoFTools::make_sparsity_pattern function you need to make 
> sure to keep the inhomogenities:
>
> DoFTools::make_sparsity_pattern(dof_handler,
>   dsp,
>   constraints,
>   /*keep_constrained_dofs = */ 
> true);
>
> The same in the assembly routine:
>
> cell->get_dof_indices(local_dof_indices);
> constraints.distribute_local_to_global(local_matrix,
> 
> local_rhs_local,
> 
> local_dof_indices,
> 

Re: [deal.II] MUMPS with PETScWrappers::MPI::BlockSparseMatrix

[Konrad]

Thank you, David, so I will see if I can maybe still use it in the Schur 
complement somehow.

Best,
Konrad

On Friday, September 6, 2019 at 7:46:23 PM UTC+2, David Wells wrote:
>
> Hi Konrad,
>
> I don't think that it is possible to use MUMPS with a block matrix for 
> exactly this reason. I think that if you want to use MUMPS you will need to 
> copy the block matrix into a non-block matrix.
>
> Thanks,
> David
>
> On Fri, Sep 6, 2019 at 12:25 PM Konrad > 
> wrote:
>
>> Dear deal.ii community,
>>
>> is it possible to use MUMPS with a PETScWrappers::MPI::BlockSparseMatrix? 
>> Don't find anything but I see that PETScWrappers::MPI::BlockSparseMatrix 
>> does not inherit from PETScWrappers::MatrixBase.
>>
>> 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 dea...@googlegroups.com .
>> To view this discussion on the web visit 
>> https://groups.google.com/d/msgid/dealii/290c3820-4c09-4515-a2f9-847691627427%40googlegroups.com
>>  
>> <https://groups.google.com/d/msgid/dealii/290c3820-4c09-4515-a2f9-847691627427%40googlegroups.com?utm_medium=email_source=footer>
>> .
>>
>

-- 
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/da7bcf06-9daa-4b52-9230-956ef11e68d3%40googlegroups.com.

[deal.II] Error with Petsc and and multithreading?

[Konrad]

Dear deal.ii community,

I am encountering an error when running the following program (only 
relevant parts):

*class NedRTStd*
*{*
*public:*
* struct Parameters*
* {*
*  // Some parameters go here*
* };*
* NedRTStd () = delete;*
* NedRTStd (Parameters );*
* ~NedRTStd ();*

* void run ();*

*private:*
* void setup_grid ();*
* void setup_system_matrix ();*
* void setup_constraints ();*
* void assemble_system ();*
* void solve_direct ();*
* void solve_iterative ();*
* void solve_iterative_preconditioned ();*
* void output_results () const;*

* MPI_Comm mpi_communicator;*

* Parameters *

* parallel::distributed::Triangulation<3> triangulation;*

* FESystem<3>fe; // Will hold 2 blocks (Nedelec-Raviart-Thomas 
pair)*
* DoFHandler<3>  dof_handler;*

* std::vector owned_partitioning;*
* std::vector relevant_partitioning;*

* AffineConstraints constraints;*

* LA::MPI::BlockSparseMatrix system_matrix;*

* LA::MPI::BlockVectorlocally_relevant_solution;*
* LA::MPI::BlockVectorsystem_rhs;*

* ConditionalOStream pcout;*
* TimerOutputcomputing_timer;*
*};*


*NedRTStd::NedRTStd (Parameters _)*
*:*
*mpi_communicator(MPI_COMM_WORLD),*
*parameters(parameters_),*
*triangulation(mpi_communicator,*
*   typename Triangulation<3>::MeshSmoothing(*
* Triangulation<3>::smoothing_on_refinement |*
* Triangulation<3>::smoothing_on_coarsening)),*
*fe (FE_Nedelec<3>(1), 1,*
* FE_RaviartThomas<3>(1), 1),*
*dof_handler (triangulation),*
*pcout(std::cout, (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)),*
*computing_timer(mpi_communicator,*
*pcout,*
*TimerOutput::summary,*
*TimerOutput::wall_times)*
*{*
*}*


*void NedRTStd::setup_system_matrix ()*
*{*
* TimerOutput::Scope t(computing_timer, "system and constraint setup");*

* dof_handler.distribute_dofs (fe);*

* if (parameters.renumber_dofs)*
* {*
* DoFRenumbering::Cuthill_McKee (dof_handler);*
* }*

* DoFRenumbering::block_wise (dof_handler);*

*//*
*// fe has 2 blocks and 6 components*
*//*
* std::vector dofs_per_block (2);*
* DoFTools::count_dofs_per_block (dof_handler, dofs_per_block);*
* const unsigned int n_sigma = dofs_per_block[0],*
*n_u = dofs_per_block[1];*

* pcout << "Number of active cells: "*
* << triangulation.n_global_active_cells()*
* << " (on " << triangulation.n_levels() << " levels)"*
* << std::endl*
* << "Number of degrees of freedom: "*
* << dof_handler.n_dofs()*
* << " (" << n_sigma << '+' << n_u << ')'*
* << std::endl;*

* owned_partitioning.resize(2);*
* owned_partitioning[0] = dof_handler.locally_owned_dofs().get_view(0, 
n_sigma);*
* owned_partitioning[1] = 
dof_handler.locally_owned_dofs().get_view(n_sigma, n_sigma + n_u);*

* IndexSet locally_relevant_dofs;*
* DoFTools::extract_locally_relevant_dofs(dof_handler, 
locally_relevant_dofs);*
* relevant_partitioning.resize(2);*
* relevant_partitioning[0] = locally_relevant_dofs.get_view(0, n_sigma);*
* relevant_partitioning[1] = locally_relevant_dofs.get_view(n_sigma, 
n_sigma + n_u);*

* setup_constraints (); // Only hanging nodes here*

* {*
* BlockDynamicSparsityPattern dsp(dofs_per_block, dofs_per_block);*

* DoFTools::make_sparsity_pattern( dof_handler, dsp, constraints, false);*
* SparsityTools::distribute_sparsity_pattern(dsp,*
* dof_handler.locally_owned_dofs_per_processor(),*
* mpi_communicator,*
* locally_relevant_dofs);*

* system_matrix.clear();*
* system_matrix.reinit(owned_partitioning, dsp, mpi_communicator);*
* }*

*//*
*// Here something happens...*

*//*
* locally_relevant_solution.reinit(owned_partitioning,*
* relevant_partitioning,*
* mpi_communicator);*

*//*
*// Here also*

*//*
* system_rhs.reinit(owned_partitioning, mpi_communicator);*
*}*

I get the following exception:

 
An error occurred in line <122> of file 
 in 
function 
   dealii::PETScWrappers::VectorBase::VectorBase() 
The violated condition was:  
   MultithreadInfo::is_running_single_threaded() 
Additional information:  
   PETSc does not support multi-threaded access, set the thread limit to 1 
in MPI_InitFinalize(). 

.
.
.


I don't have a clue what happened here since in tutorial step-55 a similar 
thing works (I can not tell what I am doing wrong). I do not invoke any 
threading deliberately.

Sorry for posting this if inappropriate - you are of course not responsible 
for debugging my code - but any hint w

Re: [deal.II] Poisson equation, implementation of boundary values, mixed finite element, semiconductor devices

[Konrad]

Hi Konrad,

I may add as a comment (without going into details): Your problem is - as 
you correctly describe above - a mixed Poisson problem. You essentially 
have two options to solve it and the way you choose influences the way you 
need to treat the boundary conditions: in primal form and in mixed form. In 
your formulation you have natural boundary conditions (Dirichlet) on Phi 
which means that you need to prescribe them in the variational form. This 
will lead to a boundary integral on the relevant part of the boundary and 
result in an additional term on the right-hand side. The boundary condition 
for E is a flux condition and in your mixed form this (Neumann) condition 
is essential. This means that you need to enforce them on the 
Raviart-Thomas space as a hard constraints. There are essentially two 
options for this:

1. If you have given a vector field that describes E on the relevant part 
of the boundary you can use the function
*VectorTools::project_boundary_values_div_conforming. *This will project 
the given vector field onto its normal component (since this is the one you 
must describe with RT elements in H(div)) and use the projected value as 
boundary condition for the normal flux of E

2. If you only are given the scalar value of the normal flux of E then you 
may have to set the relevant dofs in the RT-space by hand. This you can do 
through constraints that way:

AffineConstraints   constraints;

constraints.clear ();

DoFTools::make_hanging_node_constraints (dof_handler, constraints);

const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
const unsigned int   dofs_per_face   = fe.dofs_per_face;
std::vector local_dof_indices (dofs_per_cell);
std::vector local_dof_face_indices 
(dofs_per_face);
   
typename DoFHandler::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
cell->get_dof_indices(local_dof_indices);

for (unsigned int face_n=0;
face_n::faces_per_cell;
++face_n)
{
cell->face(face_n)->get_dof_indices(local_dof_face_indices);
if (cell->at_boundary(face_n))
{
if (cell->face(face_n)->boundary_id()==id_bottom_boundary)
{
for (unsigned int i = 0; iget_dof_indices(local_dof_indices);
constraints.distribute_local_to_global(local_matrix,

local_rhs_local,

local_dof_indices,

system_matrix,
system_rhs,
/* use 
inhomogeneities for rhs */ true);

Keep in mind that the dofs for RT elements are edge moments in 2D and 
face moments in 3D (and not nodal values).

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/48e77a20-69a8-4108-a5bc-d61adf3fba95%40googlegroups.com.

Re: [deal.II] Poisson equation, implementation of boundary values, mixed finite element, semiconductor devices

[Konrad Wiśniewski]

Wolfgang,

thank you very much once again! 
Yes, on the picture is Ex, and potential was exactly as you described. 
Today I was able to correct my calculations. I changed the right hand side 
boundary condition to be only E=0 (what is, I think, physically correct 
thing to do), and after I get the proper shape of a solution I was able to 
deduce where my mistake was. 

It looks like even the writing about a problem on forum is beneficial - it 
somehow structuring the problem. 

Thank you very much for your time! Your advises and answers really helped 
me!

Now I will start to test continuity equation so it is not excluded that I 
will ask for help in the near future on this forum (although I hope it 
won't be necessary :). 

best regards,
Konrad

W dniu czwartek, 5 września 2019 18:40:20 UTC+2 użytkownik Wolfgang 
Bangerth napisał:
>
>
> Konrad, 
>
> > I tried to give some additional constrains for electirc field E (I must 
> admit 
> > that it wasn't good idea) to avoid another unphysical results. General 
> > physical poin of view is that I should have this kind of a situation: 
> > 
> > str.png 
> > 
> > 
> > The Ex in the middle should go from 0 to -1e7 or -1e6 in the middle of 
> device 
> > and then go back to zero 
> > but my results (for Ex) look rather like this: 
> > 
> > results.png 
> > 
> > 
> > Where all Ex values, over all domain is almost -9e7 and it drops a 
> little in 
> > the middle (as it should, but not from such a big value). 
>
> What is it that you show in your figure? The x-component of the electrical 
> field? Can you also show the potential? Recall that the x-component of the 
> electrical field is something like the x-derivative of the potential Phi, 
> so 
> it seems to me like you should have a linearly decreasing potential with a 
> large gradient. Does the potential satisfy the boundary conditions you 
> impose? 
>
>
> > I couldn't find any mistake in my implementation so I've dicided to 
> change 
> > boundary condition to enforce electric filed E to be 0. 
> This is some more general advice: If you get to a situation that you don't 
> understand (electric field looks wrong) and you "fix" it by adding 
> conditions 
> you have no idea whether they are correct/physical/mathematically allowed 
> (adding boundary conditions for E), then in all likelihood you've just 
> added a 
> second bug to the first one. You really should get in the habit of 
> investigating and *understanding* what is going wrong when when something 
> looks odd, rather than *papering over* the issue. In the long run, this 
> will 
> save you time; it will also make you a better programmer and computational 
> scientists if you build a mental toolbox for how to debug problems. 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@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/07ef8e09-1c2c-435d-9ba9-d8d3e95af02b%40googlegroups.com.

[deal.II] MUMPS with PETScWrappers::MPI::BlockSparseMatrix

[Konrad]

Dear deal.ii community,

is it possible to use MUMPS with a PETScWrappers::MPI::BlockSparseMatrix? 
Don't find anything but I see that PETScWrappers::MPI::BlockSparseMatrix 
does not inherit from PETScWrappers::MatrixBase.

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/290c3820-4c09-4515-a2f9-847691627427%40googlegroups.com.

[deal.II] Function parser for tensor function

[Konrad]

Dear deal.ii community,

I was wondering if there is a simple way to parse input (e.g., from a 
parameter file) and use it similarly to the FunctionParser class (which 
works for vector valued functions) just for a TensorFunction.

Anyone did that or something similar?

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/4108d4ed-e579-4ca5-967f-3221067e7a3a%40googlegroups.com.

Re: [deal.II] Poisson equation, implementation of boundary values, mixed finite element, semiconductor devices

[Konrad Wiśniewski]

Than you for your response!

and if I may ask for some advise:

I tried to give some additional constrains for electirc field E (I must 
admit that it wasn't good idea) to avoid another unphysical results. 
General physical poin of view is that I should have this kind of a 
situation:

[image: str.png]
The Ex in the middle should go from 0 to -1e7 or -1e6 in the middle of 
device and then go back to zero 
but my results (for Ex) look rather like this:

[image: results.png]
Where all Ex values, over all domain is almost -9e7 and it drops a little 
in the middle (as it should, but not from such a big value). 

I couldn't find any mistake in my implementation so I've dicided to change 
boundary condition to enforce electric filed E to be 0. 

Do you have an intuition what I'm making wrong? 
Probably I can subtract the values on the boundary from all of the results 
but I don't think it is proper way to deal with this problem. 


W dniu środa, 4 września 2019 18:17:23 UTC+2 użytkownik Wolfgang Bangerth 
> napisał:
>>
>>
>> Konrad, 
>>
>> > I am trying to solve 2D transient problem in semiconductor devices and 
>> > I'm stuck with application of Dirichlet boundary values via 
>> > ConstrainMatrix or AffineMatrix (in the last version of this library) 
>> > 
>> > 
>> > The program (which I am rewriting) solves consecutively two equation: 
>> > (i) for a given charge density it solves Poisson equation, and than 
>> (ii) 
>> > the continuity equation to find the densities in next step. 
>> > 
>> > The poisson equation looks like this: 
>> > 
>> > 
>> > eq.png 
>> > 
>> > This is a formulation of mixed finite element problem with phi as 
>> > electric potential and E = (Ex, Ey) as a vector of electric field. 
>>
>> I suspect you have already found that out (even though you don't state 
>> it in your email) that this is *exactly* the problem step-20 solves. 
>>
>>
>> > The geometry (as well as boundary conditions) is quite simple: 
>>
>> This doesn't work. You can't impose boundary conditions on *both* Phi 
>> and E at the left and right sides. It would be equivalent to imposing 
>> both Dirichlet and Neumann values for the non-mixed Laplace equation. 
>>
>>
>> > I have choosen this kind of boundary conditions because: 
>> > 
>> > - the potential is applied on the left side of device 
>> > 
>> > - Interesting things (non-zero carrier densities) shound appeard only 
>> > near the interface of two domains  so I assumed that Ex=0, Ey=0 on 
>> > Dirichlet boundaries 
>>
>> This assumption is unphysical. From your first equation, it is clear 
>> that (up to some constant), E=-grad Phi. So if E=0 at the boundary, 
>> you'd have grad Phi=0 AND ALSO Phi=something on these boundaries. 
>>
>>
>>
>> > Unfortunately when I look at the output the E value on Dirichlet 
>> > boundary only x-component Ex is 0, but not Ey. 
>>
>> That's because for the R-T element, the degrees of freedom at a face are 
>> only the normal component of the vector. That's the only thing you can 
>> constrain, i.e., you can only prescribe E.n=something, but not the 
>> tangential component. 
>>
>> Best 
>>   W. 
>>
>>
>> -- 
>>  
>> Wolfgang Bangerth  email: bang...@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/e4256a3e-b1a3-4c90-86f0-7752dfad99fd%40googlegroups.com.

[deal.II] Re: Little question about constructor

[Konrad]

Hi Yuesu,

I think this is rather a general C++ question. Once you create an object 
all members are initialized by invoking a constructor. If you do not 
declare any constructor the compiler will silently generate a set of 
constructors. The default constructor that takes no arguments is one of 
them. This will initialize all members in the order they are given. If you 
have a member of a primitive type like double or int memory for an int or 
double will be reserved (but not initialized with a value yet). If you have 
a class member it also needs to be initialized somehow so the compiler will 
look for a default constructor.

What happens in your example is that the triangulation actually does get 
initialized by invoking its default constructor which initializes an empty 
triangulation. The dof_handler then can be called with the triangulation as 
an argument to store a (smart) pointer to the triangulation.

See this documentation for example:

https://docs.microsoft.com/en-us/cpp/cpp/constructors-cpp?view=vs-2019#constructors_in_composite_classes

It is often good to know what the C++ compiler (secretly) does without 
explicitly telling you and what functions it generates if you do not 
provide them explicitly (and when you must do so). I really recommend the 
book by Scott Meyers "Effective C++: 55 ways ..."

Hope that helps.

Best,
Konrad


On Wednesday, September 4, 2019 at 4:44:24 AM UTC+2, yuesu jin wrote:
>
> Hi all,
>
>   I have a question about the constructor initialization list. In some 
> tutorials, we initialize the dof_handler in the constructor with parameter 
> triangulation, however, the triangulation is initialized in some member 
> function within this class, which means when we initialize the class, we 
> have not had the triangulation yet, I don't know why this method works?
>
> Best regards,
> Yuesu 
>

-- 
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/fd3807d1-51e6-4daa-8eed-5da2eec1336b%40googlegroups.com.

[deal.II] Poisson equation, implementation of boundary values, mixed finite element, semiconductor devices

[Konrad Wiśniewski]

 

 Hi dear deal.ii community!


I am trying to solve 2D transient problem in semiconductor devices and I'm 
stuck with application of Dirichlet boundary values via ConstrainMatrix or 
AffineMatrix (in the last version of this library)


The program (which I am rewriting) solves consecutively two equation: (i) 
for a given charge density it solves Poisson equation, and than (ii) the 
continuity equation to find the densities in next step.

The poisson equation looks like this:

[image: eq.png]

This is a formulation of mixed finite element problem with phi as electric 
potential and E = (Ex, Ey) as a vector of electric field. 

Program uses a *Raviart-Thomas* finite element to describe E, and *Discontinous 
Galerkin* (FE_DGQ) for potential. 

Weak fe formulation:

[image: phi.png]

[image: E.png]

The geometry (as well as boundary conditions) is quite simple:

[image: bc.png]

I have choosen this kind of boundary conditions because:

- the potential is applied on the left side of device

- Interesting things (non-zero carrier densities) shound appeard only near 
the interface of two domains  so I assumed that Ex=0, Ey=0 on Dirichlet 
boundaries

- there should be no current through top and bottom of domains, so  I 
decided in most simple case to have zero electric field perpendicular to 
the boundary.

I hope those are reasonable bc. 

My implementation of boundary conditions goes like this:
For a given FESystem fe and DoFHanfler dof_handler:
dof_handler.distribute_dofs(fe);
DoFRenumbering::component_wise(dof_handler);

// make hanging node constraints
constraints.clear();
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
// add constaints to Poisson equation
const FEValuesExtractors::Vector ElectricField(0);
ComponentMask electric_field_mask = fe.component_mask(ElectricField);

DoFTools::make_zero_boundary_constraints(dof_handler,
Dirichlet, // DIRICHLET BOUNDARY INDICATOR
constraints,
electric_field_mask);

std::set no_flux_boundary;
no_flux_boundary.insert(Neumann);
VectorTools::compute_no_normal_flux_constraints(dof_handler,
0,
no_flux_boundary,
constraints);

constraints.close();

during the procedure of assembling the global matrix I use in each loop:
constraints.distribute_local_to_global( data.local_matrix,
data.local_dof_indices,
Poisson_object.system_matrix);

and for right hand side:
constraints.distribute_local_to_global(data.local_rhs,
   data.local_dof_indices,
   Poisson_object.system_rhs);

Unfortunately when I look at the output the E value on Dirichlet boundary 
only x-component Ex is 0, but not Ey. On Neumann boundary non of the 
components are 0, but of course Ey should be. Additionally the value of a 
potential (which comes from face integral on right hand side) on Dirichlet 
boundaries is not longer as it should be (eg. I ascribed 1V and in the 
output I read 1e-5V)


What did I do wrong? Is it because the Raviart-Thomas should be treated 
differently? Maybe the value that I look at in Paraview is extrapolated 
from the values from the middle of the cells face and that is why the 
values are incorrect?

Another problem is that for the simplest case all Ey on all of the domain 
should be zero. Maybe this is the problem, but I kind of don't want to 
write a code especially for the simplest case when Ey=0 -> is there a 
better way to do this?


I will appreciate any help! 

-- 
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/ff80bf4d-f269-4ad6-b131-8aab24e94f0d%40googlegroups.com.

Re: [deal.II] Distributing objects on cluster nodes according to distributed triangulation (MPI)

[Konrad]

>
> boost::mpi is a nice idea in principle, but it's not particularly well 
> designed and few people use it. Some reasons are listed here: 
>   
>
> https://scicomp.stackexchange.com/questions/3019/boostmpi-or-c-mpi-for-high-performance-scientific-applications
>  
>
> Best 
>   W. 
>
I think it is a nice interface to send not so complex objects but this 
functionality is already wrapped nicely by deal.ii. The problem starts when 
an object is not easily serializable. And, honestly, I do not have a clue 
what boost::mpi (or the boost serialization) is doing exactly when an 
object contains pointers to whatever distributed data somewhere in memory 
etc.

At any rate, Thanks for the hints.

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/0ce7a474-cc3a-4965-9c25-d310b7dc7292%40googlegroups.com.

Re: [deal.II] Distributing objects on cluster nodes according to distributed triangulation (MPI)

[Konrad]

Thank you, Wolfgang, problem solved. But in a much simpler way. My basis 
object is not easily serializable so I have to recreate it and destroy the 
old one anyway (for now) once cells are moved to other nodes. So I do not 
send send them but rather create them and work with them on each node (if 
the corresponding cell is locally owned) and just collect some serializable 
information on a master process with mpi_gather.

I have a little technical question though: Can I attach the 
mpi_communicators that you use in the deal.ii tutorials also somehow to a 
boost::mpi::communicator?

https://www.boost.org/doc/libs/1_48_0/doc/html/boost/mpi/communicator.html

Never tried it but according to the documentation it seems possible.

Best,
Konrad

On Friday, August 23, 2019 at 6:35:18 PM UTC+2, Wolfgang Bangerth wrote:
>
> On 8/23/19 7:33 AM, Konrad wrote: 
> > 
> > Thanks, I see this class makes it possible to identify cells globally in 
> a 
> > distributed triangulation. I don't see though how the objects in the 
> std::map 
> > are distributed among cluster nodes as the cells are. Am I missing 
> something here? 
>
> Yes, the data transfer is difficult because you need to figure out who 
> owns a 
> cell for which you have an entry in your local map after 
> refinement/rebalancing. Then you need to pack the data, send it to the new 
> owner, and possibly receive information about other cells you now own 
> after 
> refinement/rebalancing. 
>
> You might want to look at how other classes do this that store non-trivial 
> data on each cell. An example is ContinuousQuadratureDataTransfer. 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@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/3daf347a-3661-4eb9-ae05-6a8134141247%40googlegroups.com.

Re: [deal.II] Distributing objects on cluster nodes according to distributed triangulation (MPI)

[Konrad]

Hi Daniel,
 

> You could create an std::map using CellId (
> https://www.dealii.org/developer/doxygen/deal.II/classCellId.html) as key 
> if you need to identify cells uniquely globally.
>  
> Best,
> Daniel
>

Thanks, I see this class makes it possible to identify cells globally in a 
distributed triangulation. I don't see though how the objects in the 
std::map are distributed among cluster nodes as the cells are. Am I missing 
something here?

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/c6c93a6a-c648-48ba-a91a-61b57e8ac12c%40googlegroups.com.

[deal.II] Re: getting components of Blockvector

[Konrad]

Hi Gabriel,

I think tutorial step-22 (the part about renumbering dofs) discusses what 
you need:

You could do that when setting up your system:

std::vector block_component(dim + 1, 0);
block_component[dim] = 1;
DoFRenumbering::component_wise 
<https://www.dealii.org/current/doxygen/deal.II/namespaceDoFRenumbering.html#a75c2e3cd2d4fe86e0ec209a41a005cb0>(dof_handler,
 
block_component);

You could then access the velocity component of

 BlockVector 
<https://www.dealii.org/current/doxygen/deal.II/classBlockVector.html> 
solution;

through 

Vector 
<https://www.dealii.org/current/doxygen/deal.II/classBlockVector.html> 
velocity(solution.block(0));

Every second entry (in 2D) should then be your x-component of the velocity. ... 
at least as far as I recall (so better double check that)

Best, Konrad

Hey everyone,
>
> I have the following problem
> I have a classical velocity - pressure Blockvector "soluion" and a 
> FEsystem (initialized via
> * fe(FE_Q(degree + 1), dim, FE_Q(degree), 1)*
>
> For some testing I want to get the velocity compontents of solution in 
> seperate vectors,
> i.e. a vector u_x which contains the velocity values of solution in 
> x-direction.
>
> I am not really sure yet how to get these values.
>
> My attempt is to use FEValues and quadrature points
>
>
> u_x.reinit(n_dofs_in_one_direction);
> *FEValues fe_values(fe, 
> quadrature_formula, update_values);*
> *std::vector>  velocity_values(n_q_points);*
> *const FEValuesExtractors::Vector   velocities(0);*
>
> *for (const auto  : dof_handler.active_cell_iterators())*
> * {*
> *  fe_values.reinit(cell);*
> *  fe_values[velocities].get_function_values(solution, 
> velocity_values);*
>  
> *[...]*
>
> *}*
>
> With this procedure I get the vector of Tensors velocity_values, which 
> contains the velocity-values of "solution" in every direction at the 
> quadratue points. So now I an easily get the values of velocity_values in 
> the x-direction at the quadrature points. But how do I get the values of 
> u_x at the dofs?
>
> Does somebody have an idea how to fix this, or a better way to do handle 
> this?
>
>
> Best regards
>
> Gabriel
>

-- 
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/f4b83a6c-e494-4aae-b05d-b7ea88086013%40googlegroups.com.

[deal.II] Distributing objects on cluster nodes according to distributed triangulation (MPI)

[Konrad]

Hi deal.ii community,

I have a little conceptual question. Maybe it is dumb but I am new to MPI...

Say I have a triangulation of type 
*parallel::distributed::Triangulation* that gets distributed over N 
nodes on a cluster. Each active cell should own (or point to) an object 
that makes up a (modified non-polynomial) finite element basis. That means 
the object is not so simple in that it contains objects itself (e.g., a 
triangulation, several vectors etc).

How can I distribute a set of these basis objects (after creating the 
distributed triangulation) so that each active cell and its corresponding 
basis object are on the same cluster node? 

Is there something like a distributed STL-like container whose elements I 
can initialize on different cluster nodes (for example by querying 
*cell->is_locally_owned()*) when looping over the triangulation?

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/04fd5d49-f8e3-4414-a2d8-4370a30eb027%40googlegroups.com.

Re: [deal.II] Some problems with Eclipse IDE

[Konrad]

Dear Pawan,

setting up things by hand may be annoying. I have a github repo for this 
that you can use if you like. It works (also after I updated Ubuntu) with 
the newest version of eclipse. If you have an older version of eclipse it 
may have problems with out-of-source-builds but I strongly discourage 
in-source-builds (for the tutorials that might be ok but not for more 
complex code).

https://github.com/konsim83/Deal.ii-9.1.1_Source_code_organization

(just a suggestion that I use for my stuff, no setup of unit tests 
implemented so far but that is easy to add)

Best,
Konrad

On Monday, July 22, 2019 at 11:39:37 AM UTC+2, Pawan Kumar wrote:
>
> Dear Prof. Bangerth,
>
> Thank you very much for your reply.  Till now, the internal things between 
> cmake and eclipse is really a magic for me and it was working fine.
>
> The directory that eclipse shows for includes (in the screenshots) are the 
> installed directory of deal.II. It was the same in the earlier set up also 
> & never causes any problems.
>
> Finally, I set it up manually, as mentioned here(Setting up a project 
> using deal.II by hand) and its perfect now. 
>
> Thanks & Regards
> Pawan
>

-- 
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/01fb501a-6f45-4863-8387-32112fb38db5%40googlegroups.com.

Re: [deal.II] Parallel sparse direct solver?

[Konrad]

>
> Yes, there are such interfaces. Take a look at SparseDirectUMPFACK (no 
> PETSc required) and 
>
>
> https://dealii.org/developer/doxygen/deal.II/classPETScWrappers_1_1SparseDirectMUMPS.html
>  
>
> Best 
>   W. 
>
> Thanks! The second one I was looking for!

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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Parallel sparse direct solver?

[Konrad]

Hi deal.ii community,

I am having scalability issues with solvers for H(curl)-H(div) problems and 
H(grad)-H(curl) problems in 3D with essential boundary conditions (in a 
constraint matrix). Even for moderately sized problems (200K dofs) the 
(SparseILU preconditioned Schur-complement CG) solver needs too many 
iterations - depending on the coefficients between 500 and 1. Multigrid 
in these spaces is (as far as I know problematic). 

So I was wondering if there is a PetSc interface to MUMPS (parallel sparse 
direct solver) or if anyone has an idea for more efficient solvers for such 
problems. I am not an expert in linear solvers but their efficiency is (as 
so often) crucial here.

I would be grateful for any idea...

Cheers,
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Pure Neumann BCs in mixed form

[Konrad]

Issue resolved. I used the wrong sparsity pattern. That was just a typo in 
the code. Sorry for having bothered you with this but anyway: 

Many thanks,
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Pure Neumann BCs in mixed form

[Konrad]

uld be building a 
sparsity pattern that does not include the entries you will write into due 
to constraints on degrees of freedom such as hanging nodes or periodic 
boundary conditions. In such cases, building the sparsity pattern will 
succeed, but you will get errors such as the current one at one point or 
other when trying to write into the entries of the matrix.

It seems like I do not handle building the constraints correctly but can't 
figure out how. Sorry, that was lots of (hopefully sufficiently 
informative) code and probably it is just a minor thing.

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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Pure Neumann BCs in mixed form

[Konrad]

Hi all,

I am encountered a weird thing (bug?) and I am having trouble to understand 
it. Maybe someone has an idea:

I am trying to solve a simple elliptic problem (imagine some sort of a 
simple Laplace problem for now) with pure Neumann BCs, say -\Delta u = f, 
n\grad u = g. Now theory tells me that the volume average over f must be 
the same as the boundary average over g. If I solve that problem (with 
deal.ii) in the standard form above the Neumann condition is natural and 
enters the right hand side. I solve and I see what I expect. Fine, 
intuition confirmed.

Now I am solving it in mixed form with RT0-DGQ0 elements which are stable. 
My Neumann BC becomes essential and enters as hard constraint into the 
system but the compatibility condition should be the same. When I solve 
that system now I see things that totally do not match what I see in the 
standard form and I do not have a clue why. 

I can see from tests that also the standard form is very sensitive to 
little mismatches in the compatibility condition so I was wondering if this 
issue becomes more of a problem in mixed form.

(Btw, if I regularize the problem slightly by changing the left hand side 
-\Delta u to -\Delta u + \epsilon u  and setting f to zero (compatibility 
is not necessary then) I see something that is close to what I would expect 
in the pure Laplace case.)

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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: advection diffusion with periodic boundaries

[Konrad]

> Hi Konrad,
>
> I figured out my error. I was imposing constraints differently on the 
> advection matrix than from the other matrices. I fixed it by not using the 
> constraints.local_to_global() function and just computing the full matrix 
> and using constraints.condense() on the system matrix after adding them 
> together. 
>

Yes, that makes sense. I am glad you solved it :-)

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.
For more options, visit https://groups.google.com/d/optout.