Re: [deal.II] Using data from external matrix as source for calculations done with deal.II - best approach?

2020-02-27 Thread 'Maxi Miller' via deal.II User Group
After additional tests I noticed that I'm not only overwriting the 
target_component-component with the new values, but also the other 
components, but with 0 instead. But instead I do not want to touch them. 
How can I do that (if possible)?
Thanks!

Am Donnerstag, 27. Februar 2020 17:00:12 UTC+1 schrieb Wolfgang Bangerth:
>
>
> > | 
> > Functions::InterpolatedTensorProductGridData 
> interpolator(this->coordinate_values, 
> >   
>  this->data_storage); 
> > 
> > VectorFunctionFromScalarFunctionObject 
> interpolator_object(, 
>
> The error... 
>
> > but the compilation fails with 
> > error: ISO C++ forbids taking the address of a bound member function to 
> > form a pointer to member function.  Say 
> > ‘::Functions::InterpolatedTensorProductGridData<2>::value’. 
>
> ...comes from the line above. The problem is that 
> 
> is simply not something that C++ understands. What you need is a 
> function object that takes a Point and returns a double. You can use 
>   std::bind(::InterpolatedTensorProductGridData::value, 
> ) 
> as this argument, or maybe more modern, write 
>   [] (const Point ) { return interpolator.value(p); } 
> in place of the first constructor argument to interpolator_object. 
>
> 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/15dd3b86-8f89-412d-abec-30392e70b221%40googlegroups.com.


Re: [deal.II] Mesh refinement with periodic boundary conditions

2020-02-27 Thread Daniel Arndt
Andrew,

Did you tell the triangulation that it should take periodic boundaries into
account? My best bet (with the information given) would be that you didn't
call parallel::distributed::Triangulation::add_periodicity()

as described in https://www.dealii.org/current/doxygen/deal.II/step_45.html.

Best,
Daniel

Am Do., 27. Feb. 2020 um 17:45 Uhr schrieb Andrew Davis <
andrew@gmail.com>:

> HI all,
>
> I'm trying to implement a time-dependent convection equation using DG
> elements with periodic boundary conditions on and adaptive mesh.
>
> However, when I try to adapt the mesh using periodic using this code:
>
>   // estimate the error in each cell
>   Vector estimatedErrorPerCell(triangulation.n_active_cells());
>   //KellyErrorEstimator::estimate(dofHandler,
> QGauss(feSystem.degree + 1), {}, solution.block(0),
> estimatedErrorPerCell);
>   KellyErrorEstimator::estimate(dofHandler,
> QGauss(feSystem.degree + 1), {}, solution, estimatedErrorPerCell);
>
>   // refine cells with the largest estimated error (80 per cent of the
> error) and coarsens those cells with the smallest error (10 per cent of the
> error)
>   GridRefinement::refine_and_coarsen_fixed_fraction(triangulation,
> estimatedErrorPerCell, 0.8, 0.1);
>
>   // to prevent the decrease in time step size limit the maximal
> refinement depth of the meshlimit the maximal refinement depth of the mesh
>   if( triangulation.n_levels()>maxGridLevel ) {
> for( auto& cell :
> triangulation.active_cell_iterators_on_level(maxGridLevel) ) {
> cell->clear_refine_flag(); }
>   }
>
>   // store previous and current solution on the coarser mesh
>   std::vector > tempSolution(2);
>   tempSolution[0] = oldSolution; tempSolution[1] = solution;
>
>   // transfer the solution vectors from the old mesh to the new one
>   SolutionTransfer > transfer(dofHandler);
>
>   // prepare for the refinement---we have to prepare the transfer as well
>   *triangulation.prepare_**coarsening_and_refinement();* // for some
> reason, this doesn't tell the cells across a periodic boundary that they
> need to be refined
>   transfer.prepare_for_coarsening_and_refinement(tempSolution);
>
>   // do the refinement and reset the DoFs
>   triangulation.execute_coarsening_and_refinement();
>   SetupSystem();
>
>   std::vector > temp(2);
>   ReinitializeVector(temp[0]); ReinitializeVector(temp[1]);
>   transfer.interpolate(tempSolution, temp);
>
>   oldSolution = temp[0]; solution = temp[1];
>
>   constraints.distribute(oldSolution);
>   constraints.distribute(solution);
>
> *I get this error:*
>
> 
> An error occurred in line <992> of file
>  in function
> void dealii::{anonymous}::update_periodic_face_map_recursively(const
> typename dealii::Triangulation::cell_iterator&, const
> typename dealii::Triangulation::cell_iterator&, unsigned
> int, unsigned int, const std::bitset<3>&, std::map dealii::Triangulation::cell_iterator, unsigned int>,
> std::pair spacedim>::cell_iterator, unsigned int>, std::bitset<3> > >&) [with int dim
> = 2; int spacedim = 2; typename dealii::Triangulation spacedim>::cell_iterator = dealii::TriaIterator
> >]
> The violated condition was:
> dim == 1 || std::abs(cell_1->level() - cell_2->level()) < 2
> 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.
>
> For some reason the refinement does not seem to know that the boundaries
> are periodic and that it should refine the cells on the other side of the
> domain. I tried changing the refinement step so that instead of just
> calling triangulation.prepare_coarsening_and_refinement();, I created my
> own function that mimics the distributed triangulation (calling the
> distributed version of enforce_mesh_balance_over_periodic_boundaries)
> that does the following:
>
> while( true ) {
> triangulation.prepare_coarsening_and_refinement();
> triangulation.update_periodic_face_map();
> bool changed =
> enforce_mesh_balance_over_periodic_boundaries(triangulation);
> if( !changed ) { break; }
>   }
>
> This fixed the error but I'm getting the strange behavior---it looks like
> the cell face fluxes are no longer being 

[deal.II] Mesh refinement with periodic boundary conditions

2020-02-27 Thread Andrew Davis

HI all,

I'm trying to implement a time-dependent convection equation using DG 
elements with periodic boundary conditions on and adaptive mesh.

However, when I try to adapt the mesh using periodic using this code:

  // estimate the error in each cell
  Vector estimatedErrorPerCell(triangulation.n_active_cells());
  //KellyErrorEstimator::estimate(dofHandler, 
QGauss(feSystem.degree + 1), {}, solution.block(0), 
estimatedErrorPerCell);
  KellyErrorEstimator::estimate(dofHandler, 
QGauss(feSystem.degree + 1), {}, solution, estimatedErrorPerCell);

  // refine cells with the largest estimated error (80 per cent of the 
error) and coarsens those cells with the smallest error (10 per cent of the 
error)
  GridRefinement::refine_and_coarsen_fixed_fraction(triangulation, 
estimatedErrorPerCell, 0.8, 0.1);

  // to prevent the decrease in time step size limit the maximal refinement 
depth of the meshlimit the maximal refinement depth of the mesh
  if( triangulation.n_levels()>maxGridLevel ) {
for( auto& cell : 
triangulation.active_cell_iterators_on_level(maxGridLevel) ) { 
cell->clear_refine_flag(); }
  }

  // store previous and current solution on the coarser mesh
  std::vector > tempSolution(2);
  tempSolution[0] = oldSolution; tempSolution[1] = solution;

  // transfer the solution vectors from the old mesh to the new one
  SolutionTransfer > transfer(dofHandler);

  // prepare for the refinement---we have to prepare the transfer as well
  *triangulation.prepare_**coarsening_and_refinement();* // for some 
reason, this doesn't tell the cells across a periodic boundary that they 
need to be refined
  transfer.prepare_for_coarsening_and_refinement(tempSolution);

  // do the refinement and reset the DoFs
  triangulation.execute_coarsening_and_refinement();
  SetupSystem();

  std::vector > temp(2);
  ReinitializeVector(temp[0]); ReinitializeVector(temp[1]);
  transfer.interpolate(tempSolution, temp);

  oldSolution = temp[0]; solution = temp[1];

  constraints.distribute(oldSolution);
  constraints.distribute(solution);

*I get this error:*


An error occurred in line <992> of file 
 in function
void dealii::{anonymous}::update_periodic_face_map_recursively(const 
typename dealii::Triangulation::cell_iterator&, const 
typename dealii::Triangulation::cell_iterator&, unsigned 
int, unsigned int, const std::bitset<3>&, std::map::cell_iterator, unsigned int>, 
std::pair::cell_iterator, unsigned int>, std::bitset<3> > >&) [with int dim 
= 2; int spacedim = 2; typename dealii::Triangulation::cell_iterator = dealii::TriaIterator 
>]
The violated condition was:
dim == 1 || std::abs(cell_1->level() - cell_2->level()) < 2
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.

For some reason the refinement does not seem to know that the boundaries 
are periodic and that it should refine the cells on the other side of the 
domain. I tried changing the refinement step so that instead of just 
calling triangulation.prepare_coarsening_and_refinement();, I created my 
own function that mimics the distributed triangulation (calling the 
distributed version of enforce_mesh_balance_over_periodic_boundaries) that 
does the following:

while( true ) {
triangulation.prepare_coarsening_and_refinement();
triangulation.update_periodic_face_map();
bool changed = 
enforce_mesh_balance_over_periodic_boundaries(triangulation);
if( !changed ) { break; }
  }

This fixed the error but I'm getting the strange behavior---it looks like 
the cell face fluxes are no longer being computed correctly

Does anyone have any idea about what the issue could be? I know I left out 
a lot of code, but hopefully I included the important bits.

Thanks,
Andy
ReplyForward

 
 


-- 
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 

Re: [deal.II] Using data from external matrix as source for calculations done with deal.II - best approach?

2020-02-27 Thread 'Maxi Miller' via deal.II User Group
Yes, the lambda worked, thanks!

Am Donnerstag, 27. Februar 2020 17:00:12 UTC+1 schrieb Wolfgang Bangerth:
>
>
> > | 
> > Functions::InterpolatedTensorProductGridData 
> interpolator(this->coordinate_values, 
> >   
>  this->data_storage); 
> > 
> > VectorFunctionFromScalarFunctionObject 
> interpolator_object(, 
>
> The error... 
>
> > but the compilation fails with 
> > error: ISO C++ forbids taking the address of a bound member function to 
> > form a pointer to member function.  Say 
> > ‘::Functions::InterpolatedTensorProductGridData<2>::value’. 
>
> ...comes from the line above. The problem is that 
> 
> is simply not something that C++ understands. What you need is a 
> function object that takes a Point and returns a double. You can use 
>   std::bind(::InterpolatedTensorProductGridData::value, 
> ) 
> as this argument, or maybe more modern, write 
>   [] (const Point ) { return interpolator.value(p); } 
> in place of the first constructor argument to interpolator_object. 
>
> 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/ee9ae696-96f3-484b-bed1-fd7261f4655b%40googlegroups.com.


Re: [deal.II] Using data from external matrix as source for calculations done with deal.II - best approach?

2020-02-27 Thread Wolfgang Bangerth




|
Functions::InterpolatedTensorProductGridData 
interpolator(this->coordinate_values,
   
this->data_storage);

VectorFunctionFromScalarFunctionObject 
interpolator_object(,


The error...


but the compilation fails with
error: ISO C++ forbids taking the address of a bound member function to 
form a pointer to member function.  Say 
‘::Functions::InterpolatedTensorProductGridData<2>::value’. 


...comes from the line above. The problem is that
  
is simply not something that C++ understands. What you need is a 
function object that takes a Point and returns a double. You can use

 std::bind(::InterpolatedTensorProductGridData::value,
   )
as this argument, or maybe more modern, write
 [] (const Point ) { return interpolator.value(p); }
in place of the first constructor argument to interpolator_object.

Best
 W.
--

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

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

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/f1107e55-6c00-4bc2-8280-57f664c39246%40colostate.edu.


Re: [deal.II] Using data from external matrix as source for calculations done with deal.II - best approach?

2020-02-27 Thread 'Maxi Miller' via deal.II User Group
I implemented it the following way:
Functions::InterpolatedTensorProductGridData interpolator(this->
coordinate_values,
 
this->data_storage);

VectorFunctionFromScalarFunctionObject interpolator_object(&
interpolator.value,

target_component,

n_components_to_use);

VectorTools::internal::project_matrix_free(MappingQGeneric(1),

dof_handler,

hanging_node_constraints,

INTERPOLATION_EQUATION(fe.degree + gauss_size),

interpolator_object,

present_solution,
false,
(dim > 1 
?

 INTERPOLATION_EQUATION(fe_degree_to_use)
   : 
Quadrature(0)),
false);



but the compilation fails with 
error: ISO C++ forbids taking the address of a bound member function to 
form a pointer to member function.  Say 
‘::Functions::InterpolatedTensorProductGridData<2>::value’. 
Therefore, how can I approach that problem in a better way?
Thanks!

Am Donnerstag, 23. Januar 2020 15:39:23 UTC+1 schrieb Wolfgang Bangerth:
>
> On 1/23/20 7:29 AM, 'Maxi Miller' via deal.II User Group wrote: 
> > 
> > Why would you want to project when you can interpolate? That seems 
> > unnecessarily expensive :-) 
> > 
> > I would like to have the input data (which is interpolated) as an 
> additional 
> > component, which then is printed together with the result of the other 
> > calculations. Therefore I am not sure if there is a simpler approach 
> (after 
> > interpolation) to get the data into the component. I either can solve 
> the 
> > equation u = f for that, or use project(). Or is there another one? 
>
> What I'm saying is that you can solve for u_h=f in many different ways, 
> depending on how exactly you interpret the equality. In general, u_h is in 
> some finite element space, but f is not, so the two will not be equal 
> everywhere. Projection finds a u_h so that 
>(u_h,v_h) = (f,v_h)for all v_h 
> which is one way of interpreting the equality. Interpolation finds a u_h 
> so that 
>u_h(x_j) = f(x_j)  at all node locations x_j 
> which is another way of interpreting the equality. Neither is better than 
> the 
> other, but the second approach is *far* cheaper to compute. 
>
> 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/021304b7-81f5-40bc-91c9-a3f02832bd8a%40googlegroups.com.