Re: [deal.II] extract_dofs() extremely slow in for single-thread run?

['Maxi Miller' via deal.II User Group]

Another issue I have is for the following code piece, where I basically
have to rescale every part of the solution by a different factor (due to
using unitless variables in the calculation):
const ComponentMask surface_A_mask = fe.component_mask(surface_A);
const ComponentMask surface_B_mask = fe.component_mask(surface_B);
const ComponentMask surface_C_mask = fe.component_mask(surface_C);
const ComponentMask surface_D_mask = fe.component_mask(surface_D);

IndexSet A_value_mask = DoFTools::extract_dofs(dof_handler, surface_A_mask);
IndexSet B_value_mask = DoFTools::extract_dofs(dof_handler, surface_B_mask);
IndexSet C_value_mask = DoFTools::extract_dofs(dof_handler, surface_C_mask);
IndexSet D_value_mask = DoFTools::extract_dofs(dof_handler, surface_D_mask);

IndexSet locally_owned_elements = scaled_solution.locally_owned_elements();

A_value_mask.compress();
B_value_mask.compress();
C_value_mask.compress();
D_value_mask.compress();
locally_owned_elements.compress();

vector_t unified_test_vector(dof_handler.locally_owned_dofs(),
mpi_communicator);

for(auto index : scaled_solution.locally_owned_elements()) {

if(A_value_mask.is_element(locally_owned_elements.index_within_set(index)))
  {
  unified_test_vector(index) =
unitless_recalculator.scale_A_solution(scaled_solution(index));
 }

if(B_value_mask.is_element(locally_owned_elements.index_within_set(index)))
 {
  unified_test_vector(index) =
unitless_recalculator.scale_B_solution(scaled_solution(index));
 }

if(C_value_mask.is_element(locally_owned_elements.index_within_set(index)))
 {
  unified_test_vector(index) =
unitless_recalculator.scale_C_solution(scaled_solution(index));
 }

if(D_value_mask.is_element(locally_owned_elements.index_within_set(index)))
 {
  unified_test_vector(index) =
unitless_recalculator.scale_D_solution(scaled_solution(index));
 }
}

Again, extract_dofs takes most of the time here. Is there an alternative
approach for that, too?

Am Do., 3. Sept. 2020 um 08:56 Uhr schrieb 'develo...@googlemail.com' via
deal.II User Group :

> Yes, that worked, thanks!
> Now, when comparing, it takes 132 seconds for my initial version, and 0.41
> seconds for the version you proposed. Why is my approach so much more
> expensive?
> Thanks!
>
> Jean-Paul Pelteret schrieb am Mittwoch, 2. September 2020 um 20:34:56
> UTC+2:
>
>> I was close — it’s
>> fe.get_unit_support_points() //
>> https://www.dealii.org/developer/doxygen/deal.II/classFiniteElement.html#a5b35a290aa7dd7562911a92a13b11fee
>>
>> Moreover, defining fe_values as const disables the reinit()-call, and I
>> have to call cell->is_locally_owned(), instead of cell->locally_owned().
>>
>>
>> Coding in an email client is unforgiving, so I’m glad that you were able
>> to work past these errors.
>>
>> Best,
>> Jean-Paul
>>
>> On 02 Sep 2020, at 14:11, 'develo...@googlemail.com' via deal.II User
>> Group  wrote:
>>
>> I'm using FESystem for defining fe, and therefore I get
>> "unit_support_points is protected within this context". Is there another
>> way to get access to those points?
>> Moreover, defining fe_values as const disables the reinit()-call, and I
>> have to call cell->is_locally_owned(), instead of cell->locally_owned().
>> Thanks!
>>
>> Jean-Paul Pelteret schrieb am Dienstag, 1. September 2020 um 20:28:18
>> UTC+2:
>>
>>> Hi,
>>>
>>> 30-40% of the runtime for this operation seems a bit excessive to me.
>>> Maybe you should try using the FEValuesExtractors in conjunction with
>>> FEValues in order to get the solution at the quadrature points. If you want
>>> the solution at the support points, then you can just create a quadrature
>>> rule that has the quadrature points at the support points and initialise
>>> your FEValues object with it.
>>>
>>> So this bit of pseudo-code illustrates roughly how it could be done:
>>>
>>> //
>>> https://dealii.org/current/doxygen/deal.II/classQuadrature.html#ac002803e6ec722da5a0de102574110c9
>>> //
>>> https://dealii.org/current/doxygen/deal.II/classFiniteElement.html#ab4f6e0c83686b918fbb92716ead92313
>>> const Quadrature q_cell_support_points (fe.unit_support_points());
>>> // Assumes that your FE has unit support points; could also use the unit
>>> cell vertex positions, something else...
>>> const FEValues fe_values(fe, q_cell_support_points, update_values);
>>> const FEValuesExtractors::Scalar surface_A(0);
>>>
>>> Vector solution_vector = …; // If this is an MPI vector, then
>>> still use the same basic approach
>>>
>>> double max_value = std::numeric_limits::min();
>>> for (auto cell : dof_handler.active_cell_iterators())
>>> {
>>>   if (!cell->locally_owned()) continue;
>>>
>>>   fe_values.reinit(cell);
>>>
>>>   //
>>> https://dealii.org/current/doxygen/deal.II/classFEValuesViews_1_1Scalar.html#a7aa6b0275facea985f23a64ad690d9b1
>>>   std::vector q_point_solution_values
>>> (fe_values.n_quadrature_points);
>>>   fe_values[surface_A].get_function_values(solution_vector ,
>>> 

[deal.II] Scaling behavior of Matrix-Free test program

['Maxi Miller' via deal.II User Group]

I wrote a small test program for solving a non-linear equation using the 
RK4-solver implemented in deal.II, and assembling the right hand side using 
the matrix-free framework (code is attached). Afterwards I wanted to check 
the scaling behavior, after it should serve as a base for a larger program. 
Therefore I run several tests both on the development machine (i7-6700) 
with 8 threads and the high-performance machine (E5-2560 v4 
)
 
with 24 threads. Both machines were configured for using AVX-extensions in 
deal.II, and the program itself was compiled in release mode.

When running the program in both configurations, I compared the time it 
took for taking the first step in time:

Local machine:
MPI-ThreadsTBB-ThreadsTime (s)
1 8  170
2 4  40
4 2  20

HPC:
MPI-ThreadsTBB-ThreadsTime (s)
1 24840
2 12887
4 6  424
8 3  41
12   2  28
24   1  14

I do not fully understand that behavior: Why is the code so much slower on 
the E5 compared to the i7, except for 24 threads? Due to a different clock 
frequency, or newer structure (Broadwell vs Skylake)? Why is the transition 
from 1 MPI thread to 2 MPI threads on the i7 four times faster, but going 
from 2 MPI threads to 4 MPI threads only twice (which is expected)?
Similarly for the E5: Going from 1 thread to 2 threads does not speed up 
the code at all. Going from two to four threads halves the execution time 
(as expected), but going from four to eight results in a factor of ten. The 
steps afterwards follow the expected pattern again.

Are there any explanations for the observed behavior? And if not, what can 
I do for a deeper investigation?

Thanks!

-- 
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/01ac3d77-d9bf-4197-a9a9-1f220c0af696%40googlegroups.com.
/* -
 *
 * Copyright (C) 2009 - 2019 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, redistribute
 * it, and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * The full text of the license can be found in the file LICENSE.md at
 * the top level directory of deal.II.
 *
 * -
 *
 * Author: Wolfgang Bangerth, Texas A University, 2009, 2010
 * Timo Heister, University of Goettingen, 2009, 2010
 */
#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 
#include 

#include 
#include 
#include 

#include 
#include 
#include 

#include 
#include 
#include 

#include 

#include 
#include 

#include 
#include 

constexpr double time_step = 0.01;

constexpr unsigned int fe_degree= 3;
constexpr unsigned int n_q_points_1d= fe_degree + 2;
using Number= double;
constexpr unsigned int n_components			= 2;
constexpr unsigned int refinement_degree	= 6;

constexpr bool use_scale_function = true;
constexpr bool use_threads = use_scale_function;

namespace Step40
{
	using namespace dealii;

	using vector_t = LinearAlgebra::distributed::Vector;

	void print_vector(const vector_t _vector, const std::string ){
		std::ofstream data_out(filename);
		for(size_t i = 0; i < input_vector.size(); ++i)
			data_out << input_vector(i) << '\n';
		data_out.close();
	}

	template 
	class LaplaceOperator
	{
	public:
		static constexpr unsigned int n_quadrature_points_1d = n_points_1d;

		LaplaceOperator() : computing_times(6), local_time(0.){

		};

		void reinit(const Mapping &   mapping,
	const DoFHandler _handler,
	const AffineConstraints _node_constraints);

		void print_computing_times() const;

		void apply(const double  current_time,
   const vector_t ,
   vector_t &  dst) const;

		void 

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

['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] Using data from external matrix as source for calculations done with deal.II - best approach?

['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?

['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.

Re: [deal.II] Difference between laplace-application through a matrix and through get_function_gradients()

['Maxi Miller' via deal.II User Group]

The reason for getting wrong results was a forgotten '+'. I just inserted 
the values of the right hand side into the cell vector, but never added 
them up, unlike what I did with the matrix. Therefore I got wrong values 
for the right hand side. After this was a simple writing error, and (as far 
as I can see) nothing which can be used by others, I decided to delete the 
question. Therefore, the original question is not relevant anymore, either.
Thanks!

Am Mittwoch, 29. Januar 2020 23:19:50 UTC+1 schrieb Wolfgang Bangerth:
>
> On 1/20/20 7:59 AM, 'Maxi Miller' via deal.II User Group wrote: 
> > 
> > I wrote a short test program which should solve the diffusion equation 
> > using the time-stepping class, and implemented both methods. When 
> > calculating the matrix and applying it to my solution vector, I get a 
> > different result compared to reading the gradients from the solution 
> > with get_function_gradients() and multiplying it with the gradients 
> > returned by shape_grad(). The results obtained from the matrix 
> > multiplication are correct compared to the expected solution, the 
> > results obtained from the direct approach are not. Why is that? 
>
> Maxi -- if I understand you correctly, you're asking what the difference 
> is between computing 
>
>F_i = \int \grad\varphi_i  .  grad u_h 
>
> and 
>
>F_i = (AU)_i 
>
> where A=Laplace matrix, U=coefficient vector corresponding to u_h. 
>
> There shouldn't be a difference in principle, but you have to pay 
> attention to what hanging nodes and Dirichlet boundary conditions do. In 
> particular, you might have to call F.condense() in the first case. 
>
> You only say that the results are different, but not *how* they are 
> different. Have you looked at that? Are these two vectors different only 
> in hanging nodes? Only for shape functions at the boundary? 
>
> 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/327036e0-7343-488e-949e-8ce1bb42f829%40googlegroups.com.

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

['Maxi Miller' via deal.II User Group]

Am Mittwoch, 22. Januar 2020 19:23:52 UTC+1 schrieb Wolfgang Bangerth:
>
> On 1/22/20 11:17 AM, 'Maxi Miller' via deal.II User Group wrote: 
> > As source term for calculations I am using data which is written in an 
> > external matrix (because it has been calculated by an external 
> > function), together with the x- and y-values for each matrix entry. My 
> > current approach for using this data is to use an interpolation function 
> > from GSL, and call it for each cell while looping over them. This gives 
> > me the value at that position, interpolated. 
>
> You might also want to look at the InterpolatedTensorProductGridData and 
> InterpolatedUniformGridData classes. Since these are derived from 
> Function, you can send them into VectorTools::project if you want. 
>
>
> > Now, though, I would like to improve that approach. I was thinking about 
> > using the project-function and interpolating the matrix onto a component 
> > in my FESystem, but as far as I can see the project-function does not 
> > allow me to only change one component in my solution while leaving the 
> > other components untouched. Of course, I also could loop once over all 
> > cells before starting the main loop and solve the equation u = f, with f 
> > the value at each position, but I am not sure how efficient that is, 
> > especially when I can use the matrix-free vectorized approach for 
> > solving for the other components. 
>
> 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? 

>
> But regardless, if you can't find a way to deal with just component, 
> embed your data into a vector-valued function where the other components 
> are all zero, and then project/interpolate the whole vector-valued 
> function onto your FE space. You might be interested in the 
> VectorFunctionFromScalarFunctionObject class for this purpose. 
>
> 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/1156c2ae-4180-4446-9a5f-e167a01e9c44%40googlegroups.com.

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

['Maxi Miller' via deal.II User Group]

As source term for calculations I am using data which is written in an 
external matrix (because it has been calculated by an external function), 
together with the x- and y-values for each matrix entry. My current 
approach for using this data is to use an interpolation function from GSL, 
and call it for each cell while looping over them. This gives me the value 
at that position, interpolated.
Now, though, I would like to improve that approach. I was thinking about 
using the project-function and interpolating the matrix onto a component in 
my FESystem, but as far as I can see the project-function does not allow me 
to only change one component in my solution while leaving the other 
components untouched. Of course, I also could loop once over all cells 
before starting the main loop and solve the equation u = f, with f the 
value at each position, but I am not sure how efficient that is, especially 
when I can use the matrix-free vectorized approach for solving for the 
other components.
Therefore, are there other (maybe more efficient) approaches which I can 
use for transferring and interpolating the values from the original matrix 
to something which I then can use in my deal.II-functions?
Thanks!

-- 
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/290765bf-ae24-404d-9e36-64cb887a1e22%40googlegroups.com.

[deal.II] Usage of vdSin/similar functions from the MKL together with VectorizedArray

['Maxi Miller' via deal.II User Group]

As described in step-48 it is possible to use vectorized functions, such as 
the function vdSin() from the MKL-library from Intel, together with 
VectorizedArray-variables. Unfortunately I could not find any example about 
how I can do that. Is the approach to simply use a code similar to 
VectorizedArray a, b;
//fill a;
std::vector a_tmp(VectorizedArray::n_array_elements), b_tmp(
VectorizedArray::n_array_elements);
a.store(a_tmp.data());

vdSin(VectorizedArray::n_array_elements, _tmp, _tmp);
b.load(b_tmp.data());

or is there a simpler approach possible?

-- 
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/8678f432-7234-45eb-b508-77d3792402b3%40googlegroups.com.

[deal.II] Re: Implementation of a solver for the non-linear diffusion equation using the TimeStepping-Class

['Maxi Miller' via deal.II User Group]

Found my mistake, I was iterating over the wrong things.

Am Sonntag, 19. Januar 2020 14:54:16 UTC+1 schrieb Maxi Miller:
>
> I added a small MWE-code to the question, to make debugging easier. 
> Switching between the linear heat equation and the non-linear heat equation 
> can be done using the constexpr-variable use_nonlinear in line 79.
> Furthermore, I changed the solution from exp(-pi^2*t) to exp(-2 * pi^2 * 
> t), and changed the value of f accordingly.
> Thanks!
>
> Am Samstag, 18. Januar 2020 12:11:43 UTC+1 schrieb Maxi Miller:
>>
>> I tried to implement a solver for the non-linear diffusion equation 
>> (\partial_t u = grad(u(grad u)) - f) using the TimeStepping-Class, the 
>> EmbeddedExplicitRungeKutta-Method and (for assembly) the matrix-free 
>> approach. For initial tests I used the linear heat equation with the 
>> solution u = exp(-pi * pi * t) * sin(pi * x) * sin(pi * y) and the assembly 
>> routine
>> template 
>> void LaplaceOperator::local_apply_cell(
>> const MatrixFree &   data,
>> vector_t &  dst,
>> const vector_t ,
>> const std::pair & cell_range) 
>> const
>> {
>> FEEvaluation> Number> phi(data);
>>
>> for (unsigned int cell = cell_range.first; cell < cell_range.
>> second; ++cell)
>> {
>> phi.reinit(cell);
>> phi.read_dof_values_plain(src);
>> phi.evaluate(false, true);
>> for (unsigned int q = 0; q < phi.n_q_points; ++q)
>> {
>> auto gradient_coefficient = 
>> calculate_gradient_coefficient(phi.get_gradient(q));
>> phi.submit_gradient(gradient_coefficient, q);
>> }
>> phi.integrate_scatter(false, true, dst);
>> }
>> }
>>
>>
>> and 
>> template 
>> inline DEAL_II_ALWAYS_INLINE
>> Tensor<1, n_components_to_use, Tensor<1, dim, Number>> 
>> calculate_gradient_coefficient(
>> #if defined(USE_NONLINEAR) || defined(USE_ADVECTION)
>> const Tensor<1, n_components_to_use, Number> _value,
>> #endif
>> const Tensor<1, n_components_to_use, Tensor<1, dim, Number>> 
>> _gradient){
>> Tensor<1, n_components_to_use, Tensor<1, dim, Number>> ret_val;
>> for(size_t component = 0; component < n_components_to_use; ++
>> component){
>> for(size_t d = 0; d < dim; ++d){
>> ret_val[component][d] = -1. / (2 * M_PI * M_PI) * 
>> input_gradient[component][d];
>> }
>> }
>> return ret_val;
>> }
>>
>>
>> This approach works, and delivers correct results. Now I wanted to test 
>> the same approach for the non-linear diffusion equation with f = -exp(-2 * 
>> pi^2 * t) * 0.5 * pi^2 * (-cos(2 * pi * (x - y)) - cos(2 * pi * (x + y)) + 
>> cos(2 * pi * x) + cos(2 * pi * y)), which should be the solution to grad(u 
>> (grad u)) with u = exp(-pi^2*t) * sin(pi * x) * sin(pi * y). Thus, I 
>> changed the routines to
>> template 
>> void LaplaceOperator::local_apply_cell(
>> const MatrixFree &   data,
>> vector_t &  dst,
>> const vector_t ,
>> const std::pair & cell_range) 
>> const
>> {
>> FEEvaluation> Number> phi(data);
>>
>> for (unsigned int cell = cell_range.first; cell < cell_range.
>> second; ++cell)
>> {
>> phi.reinit(cell);
>> phi.read_dof_values_plain(src);
>> phi.evaluate(true, true);
>> for(size_t q = 0; q < phi.n_q_points; ++q){
>> auto value = phi.get_value(q);
>> auto gradient = phi.get_gradient(q);
>> phi.submit_value(calculate_value_coefficient(value,
>>  phi.
>> quadrature_point(q),
>>  local_time), 
>> q);
>> phi.submit_gradient(calculate_gradient_coefficient(value,
>>   
>>  gradient), q);
>> }
>> phi.integrate_scatter(true, true, dst);
>> }
>> }
>>
>> and 
>> template 
>> inline DEAL_II_ALWAYS_INLINE
>> Tensor<1, n_components_to_use, Number> calculate_value_coefficient(
>> const Tensor<1, n_components_to_use, Number> _value,
>>
>> const Point ,
>>
>> const double ){
>> Tensor<1, n_components_to_use, Number> ret_val;
>> //(void) input_value;
>> (void) input_value;
>> for(size_t component = 0; component < n_components_to_use; ++
>> component){
>> const double x = point[component][0];
>> const double y = point[component][1];
>> ret_val[component] = (- exp(-2 * M_PI * M_PI * time)
>>

Re: [deal.II] Re: Application of inverse mass matrix in matrix-free context using cell_loop instead of scale()

['Maxi Miller' via deal.II User Group]

Found a bug in the library, reported here: 
https://github.com/dealii/dealii/issues/9405, and therefore I got wrong 
results.

Am Sonntag, 19. Januar 2020 13:50:58 UTC+1 schrieb Maxi Miller:
>
> Hei,
> I attached a working MWE. Changing between the approach suggested above 
> and scale() is done by changing the variable use_scale_function. The output 
> contains both the expected solution and the obtained solution.
> My approach for replacing scale() with a local loop is
> data.cell_loop(::local_apply_cell,
>this,
>dst,
>src,
>[&](const unsigned int start_range, const 
> unsigned int end_range){
> for(size_t i = start_range; i < end_range; ++i){
> dst.local_element(i) = 0;
> }
> },
> [&](const unsigned int start_range, const unsigned int 
> end_range){
> for(unsigned int i = start_range; i < end_range; ++i){
> dst.local_element(i) = dst.local_element(i) * 
> inv_mass_matrix.local_element(i);
> }
> })
>
> I expect my result to follow the solution function exp(-2 * pi^2 * t) * 
> sin(pi * x) * sin(pi * y). While that works for the scale()-approach, the 
> result for integrated approach does not change from time step to time step. 
> Currently I am only running it with a single MPI thread, after encountering 
> issues when running with several threads.
> Thanks!
>
> Am Samstag, 18. Januar 2020 18:42:04 UTC+1 schrieb Daniel Arndt:
>>
>> Maxi,
>>
>> As usual, it is much easier to help if you provide a complete minimal 
>> example and say how the result differs from what you expect.
>> Does it only scale certain vector entries? Are the results correct when 
>> running with one MPI process?
>> How does your approach differ from 
>> https://github.com/dealii/dealii/blob/b84270a1d4099292be5b3d43c2ea65f3ee005919/tests/matrix_free/pre_and_post_loops_01.cc#L100-L121
>> ?
>>
>> Best,
>> Daniel
>>
>> Am Sa., 18. Jan. 2020 um 12:05 Uhr schrieb 'Maxi Miller' via deal.II User 
>> Group :
>>
>>> I tried implementing it as 
>>> data.cell_loop(::local_apply_cell,
>>>this,
>>>dst,
>>>src,
>>>//true,
>>>[&](const unsigned int start_range, const 
>>> unsigned int end_range){
>>> for(size_t i = start_range; i < end_range; ++i){
>>> dst.local_element(i) = 0;
>>> }
>>> },
>>> [&](const unsigned int start_range, const unsigned int end_range
>>> ){
>>> for(unsigned int i = start_range; i < end_range; ++i){
>>> dst.local_element(i) = dst.local_element(i) * 
>>> inv_mass_matrix.local_element(i);
>>> }
>>> });
>>>
>>> but the result was not correct. Thus, I assume I have to do something 
>>> else?
>>> Thanks!
>>>
>>> Am Samstag, 18. Januar 2020 17:12:34 UTC+1 schrieb peterrum:
>>>>
>>>> Yes, like here 
>>>> https://github.com/dealii/dealii/blob/b84270a1d4099292be5b3d43c2ea65f3ee005919/tests/matrix_free/pre_and_post_loops_01.cc#L100-L121
>>>>
>>>> On Saturday, 18 January 2020 12:57:24 UTC+1, Maxi Miller wrote:
>>>>>
>>>>> In step-48 the inverse mass matrix is applied by moving the inverse 
>>>>> data into a vector and applying the function scale(), i.e. as in the 
>>>>> following code:
>>>>> data.cell_loop(::local_apply_cell,
>>>>>this,
>>>>>dst,
>>>>>src,
>>>>>true);
>>>>> computing_times[0] += timer.wall_time();
>>>>> timer.restart();
>>>>>
>>>>>dst.scale(inv_mass_matrix);
>>>>>
>>>>>
>>>>>
>>>>> Now I would like to do the same, but use a cell_loop instead of 
>>>>> scale(). Thus, I created an additional function called 
>>>>> "local_apply_inverse_mass_matrix" as 
>>>>> template 
>>>>> void LaplaceOperator::
>>>>> local_apply_inverse_mass_matrix(
>>>>>   

[deal.II] Difference between laplace-application through a matrix and through get_function_gradients()

['Maxi Miller' via deal.II User Group]

When applying the laplace operator to my solution vector, what is the 
difference between forming an explicit matrix containing the laplace 
operator by using 
template 
void LaplaceProblem::assemble_mass_matrix (){
TimerOutput::Scope t(computing_timer, "Mass matrix assembly");
mass_matrix = 0.;
system_matrix = 0.;

const 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.dofs_per_cell;
const unsigned int n_q_points= quadrature_formula.size();
FullMatrix  cell_mass_matrix(dofs_per_cell,
 dofs_per_cell),
cell_matrix(dofs_per_cell,
dofs_per_cell);
std::vector local_dof_indices(dofs_per_cell
);

for (const auto  : dof_handler.active_cell_iterators())
if (cell->is_locally_owned())
{
cell_mass_matrix = 0.;
cell_matrix = 0.;
fe_values.reinit(cell);

for (unsigned int q = 0; q < n_q_points; ++q)
{
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
for (unsigned int j = 0; j < dofs_per_cell; ++j){
cell_matrix(i, j) += -1.
* fe_values.shape_grad (i, q)
* fe_values.shape_grad (j, q)
* fe_values.JxW (q);
cell_mass_matrix(i, j) += (fe_values.shape_value
(i, q)
   * fe_values.
shape_value(j, q))
* fe_values.JxW(q);
}
}
}
cell->get_dof_indices(local_dof_indices);

constraints.distribute_local_to_global(cell_mass_matrix,
   local_dof_indices,
   mass_matrix);
constraints.distribute_local_to_global(cell_matrix,
   local_dof_indices,
   system_matrix);
}
mass_matrix.compress(VectorOperation::add);
system_matrix.compress(VectorOperation::add);
}



and assembling the result vector directly by getting the gradients from the 
solution vector via get_function_gradients() and multiplying it afterwards 
with shape_grad(), by using
template 
void LaplaceProblem::assemble_rhs(const vector_t , const double 
local_time, vector_t ){
TimerOutput::Scope t(computing_timer, "RHS assembly");

dst.reinit(locally_owned_dofs, mpi_communicator);
dst = 0.;

const 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.dofs_per_cell;
const unsigned int n_q_points= quadrature_formula.size();
Vector cell_rhs(dofs_per_cell);
std::vector local_dof_indices(dofs_per_cell
);
std::vector> solution_gradients(n_q_points);

for (const auto  : dof_handler.active_cell_iterators())
if (cell->is_locally_owned())
{
cell_rhs = 0.;
fe_values.reinit(cell);

fe_values.get_function_gradients(y, solution_gradients);

for (unsigned int q = 0; q < n_q_points; ++q)
{
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
cell_rhs(i) = (-1.
   * solution_gradients[q]
   * fe_values.shape_grad(i, q))
* fe_values.JxW(q);
}
}
cell->get_dof_indices(local_dof_indices);

constraints.distribute_local_to_global(cell_rhs,
   local_dof_indices,
   dst);
}
dst.compress(VectorOperation::add);
}


I wrote a short test program which should solve the diffusion equation 
using the time-stepping class, and implemented both methods. When 
calculating the matrix and applying it to my solution vector, I get a 
different result compared to reading the gradients from the solution with 
get_function_gradients() and 

[deal.II] Extension of gather_evaluate for read_dof_values_plain()?

['Maxi Miller' via deal.II User Group]

Currently gather_evaluate is a replacement for the combination of 
read_dof_values() and evaluate(), but it can not be used if one has to use 
read_dof_values_plain() instead. Therefore I was wondering if it was not 
implemented simply due to a lack of time, or because it is more difficult 
to implement it for that use case?
Thanks!

-- 
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/349d80e6-9d9b-461e-9635-e9173073ec5d%40googlegroups.com.

[deal.II] Re: Implementation of a solver for the non-linear diffusion equation using the TimeStepping-Class

['Maxi Miller' via deal.II User Group]

I added a small MWE-code to the question, to make debugging easier. 
Switching between the linear heat equation and the non-linear heat equation 
can be done using the constexpr-variable use_nonlinear in line 79.
Furthermore, I changed the solution from exp(-pi^2*t) to exp(-2 * pi^2 * 
t), and changed the value of f accordingly.
Thanks!

Am Samstag, 18. Januar 2020 12:11:43 UTC+1 schrieb Maxi Miller:
>
> I tried to implement a solver for the non-linear diffusion equation 
> (\partial_t u = grad(u(grad u)) - f) using the TimeStepping-Class, the 
> EmbeddedExplicitRungeKutta-Method and (for assembly) the matrix-free 
> approach. For initial tests I used the linear heat equation with the 
> solution u = exp(-pi * pi * t) * sin(pi * x) * sin(pi * y) and the assembly 
> routine
> template 
> void LaplaceOperator::local_apply_cell(
> const MatrixFree &   data,
> vector_t &  dst,
> const vector_t ,
> const std::pair & cell_range) 
> const
> {
> FEEvaluation > phi(data);
>
> for (unsigned int cell = cell_range.first; cell < cell_range.
> second; ++cell)
> {
> phi.reinit(cell);
> phi.read_dof_values_plain(src);
> phi.evaluate(false, true);
> for (unsigned int q = 0; q < phi.n_q_points; ++q)
> {
> auto gradient_coefficient = calculate_gradient_coefficient
> (phi.get_gradient(q));
> phi.submit_gradient(gradient_coefficient, q);
> }
> phi.integrate_scatter(false, true, dst);
> }
> }
>
>
> and 
> template 
> inline DEAL_II_ALWAYS_INLINE
> Tensor<1, n_components_to_use, Tensor<1, dim, Number>> 
> calculate_gradient_coefficient(
> #if defined(USE_NONLINEAR) || defined(USE_ADVECTION)
> const Tensor<1, n_components_to_use, Number> _value,
> #endif
> const Tensor<1, n_components_to_use, Tensor<1, dim, Number>> &
> input_gradient){
> Tensor<1, n_components_to_use, Tensor<1, dim, Number>> ret_val;
> for(size_t component = 0; component < n_components_to_use; ++
> component){
> for(size_t d = 0; d < dim; ++d){
> ret_val[component][d] = -1. / (2 * M_PI * M_PI) * 
> input_gradient[component][d];
> }
> }
> return ret_val;
> }
>
>
> This approach works, and delivers correct results. Now I wanted to test 
> the same approach for the non-linear diffusion equation with f = -exp(-2 * 
> pi^2 * t) * 0.5 * pi^2 * (-cos(2 * pi * (x - y)) - cos(2 * pi * (x + y)) + 
> cos(2 * pi * x) + cos(2 * pi * y)), which should be the solution to grad(u 
> (grad u)) with u = exp(-pi^2*t) * sin(pi * x) * sin(pi * y). Thus, I 
> changed the routines to
> template 
> void LaplaceOperator::local_apply_cell(
> const MatrixFree &   data,
> vector_t &  dst,
> const vector_t ,
> const std::pair & cell_range) 
> const
> {
> FEEvaluation > phi(data);
>
> for (unsigned int cell = cell_range.first; cell < cell_range.
> second; ++cell)
> {
> phi.reinit(cell);
> phi.read_dof_values_plain(src);
> phi.evaluate(true, true);
> for(size_t q = 0; q < phi.n_q_points; ++q){
> auto value = phi.get_value(q);
> auto gradient = phi.get_gradient(q);
> phi.submit_value(calculate_value_coefficient(value,
>  phi.
> quadrature_point(q),
>  local_time), 
> q);
> phi.submit_gradient(calculate_gradient_coefficient(value,
>gradient
> ), q);
> }
> phi.integrate_scatter(true, true, dst);
> }
> }
>
> and 
> template 
> inline DEAL_II_ALWAYS_INLINE
> Tensor<1, n_components_to_use, Number> calculate_value_coefficient(
> const Tensor<1, n_components_to_use, Number> _value,
>
> const Point ,
>
> const double ){
> Tensor<1, n_components_to_use, Number> ret_val;
> //(void) input_value;
> (void) input_value;
> for(size_t component = 0; component < n_components_to_use; ++
> component){
> const double x = point[component][0];
> const double y = point[component][1];
> ret_val[component] = (- exp(-2 * M_PI * M_PI * time)
>   * 0.5 * M_PI * M_PI * (-cos(2 * M_PI * (x 
> - y))
>  - cos(2 * M_PI * 
> (x + y))
>  + cos(2 * M_PI * 
> x)
>  

Re: [deal.II] Re: Application of inverse mass matrix in matrix-free context using cell_loop instead of scale()

['Maxi Miller' via deal.II User Group]

Hei,
I attached a working MWE. Changing between the approach suggested above and 
scale() is done by changing the variable use_scale_function. The output 
contains both the expected solution and the obtained solution.
My approach for replacing scale() with a local loop is
data.cell_loop(::local_apply_cell,
   this,
   dst,
   src,
   [&](const unsigned int start_range, const 
unsigned int end_range){
for(size_t i = start_range; i < end_range; ++i){
dst.local_element(i) = 0;
}
},
[&](const unsigned int start_range, const unsigned int end_range
){
for(unsigned int i = start_range; i < end_range; ++i){
dst.local_element(i) = dst.local_element(i) * 
inv_mass_matrix.local_element(i);
}
})

I expect my result to follow the solution function exp(-2 * pi^2 * t) * 
sin(pi * x) * sin(pi * y). While that works for the scale()-approach, the 
result for integrated approach does not change from time step to time step. 
Currently I am only running it with a single MPI thread, after encountering 
issues when running with several threads.
Thanks!

Am Samstag, 18. Januar 2020 18:42:04 UTC+1 schrieb Daniel Arndt:
>
> Maxi,
>
> As usual, it is much easier to help if you provide a complete minimal 
> example and say how the result differs from what you expect.
> Does it only scale certain vector entries? Are the results correct when 
> running with one MPI process?
> How does your approach differ from 
> https://github.com/dealii/dealii/blob/b84270a1d4099292be5b3d43c2ea65f3ee005919/tests/matrix_free/pre_and_post_loops_01.cc#L100-L121
> ?
>
> Best,
> Daniel
>
> Am Sa., 18. Jan. 2020 um 12:05 Uhr schrieb 'Maxi Miller' via deal.II User 
> Group >:
>
>> I tried implementing it as 
>> data.cell_loop(::local_apply_cell,
>>this,
>>dst,
>>src,
>>//true,
>>[&](const unsigned int start_range, const unsigned 
>> int end_range){
>> for(size_t i = start_range; i < end_range; ++i){
>> dst.local_element(i) = 0;
>> }
>> },
>> [&](const unsigned int start_range, const unsigned int end_range
>> ){
>> for(unsigned int i = start_range; i < end_range; ++i){
>> dst.local_element(i) = dst.local_element(i) * 
>> inv_mass_matrix.local_element(i);
>> }
>> });
>>
>> but the result was not correct. Thus, I assume I have to do something 
>> else?
>> Thanks!
>>
>> Am Samstag, 18. Januar 2020 17:12:34 UTC+1 schrieb peterrum:
>>>
>>> Yes, like here 
>>> https://github.com/dealii/dealii/blob/b84270a1d4099292be5b3d43c2ea65f3ee005919/tests/matrix_free/pre_and_post_loops_01.cc#L100-L121
>>>
>>> On Saturday, 18 January 2020 12:57:24 UTC+1, Maxi Miller wrote:
>>>>
>>>> In step-48 the inverse mass matrix is applied by moving the inverse 
>>>> data into a vector and applying the function scale(), i.e. as in the 
>>>> following code:
>>>> data.cell_loop(::local_apply_cell,
>>>>this,
>>>>dst,
>>>>src,
>>>>true);
>>>> computing_times[0] += timer.wall_time();
>>>> timer.restart();
>>>>
>>>>dst.scale(inv_mass_matrix);
>>>>
>>>>
>>>>
>>>> Now I would like to do the same, but use a cell_loop instead of 
>>>> scale(). Thus, I created an additional function called 
>>>> "local_apply_inverse_mass_matrix" as 
>>>> template 
>>>> void LaplaceOperator::
>>>> local_apply_inverse_mass_matrix(
>>>> const MatrixFree &   data,
>>>> LinearAlgebra::distributed::Vector &  dst,
>>>> const LinearAlgebra::distributed::Vector ,
>>>> const std::pair & 
>>>> cell_range) const
>>>> {
>>>> (void) data;
>>>> (void) cell_range;
>>>> dst = src;
>>>> dst.scale(inv_mass_matrix);
>>>> }
>>>>
>>>> When running that code, though, using
>>>> LinearAlgebra::

[deal.II] Re: Application of inverse mass matrix in matrix-free context using cell_loop instead of scale()

['Maxi Miller' via deal.II User Group]

I tried implementing it as 
data.cell_loop(::local_apply_cell,
   this,
   dst,
   src,
   //true,
   [&](const unsigned int start_range, const unsigned 
int end_range){
for(size_t i = start_range; i < end_range; ++i){
dst.local_element(i) = 0;
}
},
[&](const unsigned int start_range, const unsigned int end_range){
for(unsigned int i = start_range; i < end_range; ++i){
dst.local_element(i) = dst.local_element(i) * 
inv_mass_matrix.local_element(i);
}
});

but the result was not correct. Thus, I assume I have to do something else?
Thanks!

Am Samstag, 18. Januar 2020 17:12:34 UTC+1 schrieb peterrum:
>
> Yes, like here 
> https://github.com/dealii/dealii/blob/b84270a1d4099292be5b3d43c2ea65f3ee005919/tests/matrix_free/pre_and_post_loops_01.cc#L100-L121
>
> On Saturday, 18 January 2020 12:57:24 UTC+1, Maxi Miller wrote:
>>
>> In step-48 the inverse mass matrix is applied by moving the inverse data 
>> into a vector and applying the function scale(), i.e. as in the following 
>> code:
>> data.cell_loop(::local_apply_cell,
>>this,
>>dst,
>>src,
>>true);
>> computing_times[0] += timer.wall_time();
>> timer.restart();
>>
>>dst.scale(inv_mass_matrix);
>>
>>
>>
>> Now I would like to do the same, but use a cell_loop instead of scale(). 
>> Thus, I created an additional function called 
>> "local_apply_inverse_mass_matrix" as 
>> template 
>> void LaplaceOperator::
>> local_apply_inverse_mass_matrix(
>> const MatrixFree &   data,
>> LinearAlgebra::distributed::Vector &  dst,
>> const LinearAlgebra::distributed::Vector ,
>> const std::pair & cell_range) 
>> const
>> {
>> (void) data;
>> (void) cell_range;
>> dst = src;
>> dst.scale(inv_mass_matrix);
>> }
>>
>> When running that code, though, using
>> LinearAlgebra::distributed::Vector tmp(src);
>>
>> data.initialize_dof_vector(tmp);
>> data.initialize_dof_vector(dst);
>> data.cell_loop(::local_apply_cell,
>>this,
>>tmp,
>>src,
>>true);
>> computing_times[0] += timer.wall_time();
>> timer.restart();
>>
>> data.cell_loop(::local_apply_inverse_mass_matrix,
>>this,
>>dst,
>>tmp,
>>true);
>> computing_times[1] += timer.wall_time();
>>
>> computing_times[3] += 1.;
>>
>>
>> I get the error 
>> An error occurred in line <3338> of file > matrix_free/matrix_free.h> in function
>> void dealii::internal::VectorDataExchange> VectorizedArrayType>::compress_start(unsigned int, VectorType&) [with 
>> VectorType = dealii::LinearAlgebra::distributed::Vector; typename 
>> std::enable_if<(dealii::internal::has_compress_start::value 
>> && dealii::internal::has_exchange_on_subset::value), 
>> VectorType>::type*  = 0; int dim = 2; Number = double; 
>> VectorizedArrayType = dealii::VectorizedArray]
>> The violated condition was: 
>> vec.has_ghost_elements() == false
>> 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).
>>
>> Stacktrace:
>> ---
>> #0  ./MF_FES_RK4-Test: void dealii::internal::VectorDataExchange<2, 
>> double, dealii::VectorizedArray 
>> >::compress_start> dealii::MemorySpace::Host>, 
>> (dealii::LinearAlgebra::distributed::Vector> dealii::MemorySpace::Host>*)0>(unsigned int, 
>> dealii::LinearAlgebra::distributed::Vector> dealii::MemorySpace::Host>&)
>> #1  ./MF_FES_RK4-Test: void dealii::internal::compress_start<2, 
>> dealii::LinearAlgebra::distributed::Vector> dealii::MemorySpace::Host>, double, dealii::VectorizedArray, 
>> (dealii::LinearAlgebra::distributed::Vector> dealii::MemorySpace::Host>*)0>(dealii::LinearAlgebra::distributed::Vector>  
>> dealii::MemorySpace::Host>&, dealii::internal::VectorDataExchange<2, 
>> double, dealii::VectorizedArray >&, unsigned int)
>> #2  ./MF_FES_RK4-Test: dealii::internal::MFWorker> double, dealii::VectorizedArray >, 
>> dealii::LinearAlgebra::distributed::Vector> 

[deal.II] Re: Application of inverse mass matrix in matrix-free context using cell_loop instead of scale()

['Maxi Miller' via deal.II User Group]

I.e. I should add an elementwise multiplication with the inverse mass matrix 
vector as postprocessing function?

-- 
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/200f6ca1-3a47-4e2e-bc7e-544f8d112697%40googlegroups.com.

[deal.II] Application of inverse mass matrix in matrix-free context using cell_loop instead of scale()

['Maxi Miller' via deal.II User Group]

In step-48 the inverse mass matrix is applied by moving the inverse data 
into a vector and applying the function scale(), i.e. as in the following 
code:
data.cell_loop(::local_apply_cell,
   this,
   dst,
   src,
   true);
computing_times[0] += timer.wall_time();
timer.restart();

   dst.scale(inv_mass_matrix);



Now I would like to do the same, but use a cell_loop instead of scale(). 
Thus, I created an additional function called 
"local_apply_inverse_mass_matrix" as 
template 
void LaplaceOperator::
local_apply_inverse_mass_matrix(
const MatrixFree &   data,
LinearAlgebra::distributed::Vector &  dst,
const LinearAlgebra::distributed::Vector ,
const std::pair & cell_range) 
const
{
(void) data;
(void) cell_range;
dst = src;
dst.scale(inv_mass_matrix);
}

When running that code, though, using
LinearAlgebra::distributed::Vector tmp(src);

data.initialize_dof_vector(tmp);
data.initialize_dof_vector(dst);
data.cell_loop(::local_apply_cell,
   this,
   tmp,
   src,
   true);
computing_times[0] += timer.wall_time();
timer.restart();

data.cell_loop(::local_apply_inverse_mass_matrix,
   this,
   dst,
   tmp,
   true);
computing_times[1] += timer.wall_time();

computing_times[3] += 1.;


I get the error 
An error occurred in line <3338> of file  in function
void dealii::internal::VectorDataExchange::compress_start(unsigned int, VectorType&) [with 
VectorType = dealii::LinearAlgebra::distributed::Vector; typename 
std::enable_if<(dealii::internal::has_compress_start::value && 
dealii::internal::has_exchange_on_subset::value), VectorType>::
type*  = 0; int dim = 2; Number = double; VectorizedArrayType = 
dealii::VectorizedArray]
The violated condition was: 
vec.has_ghost_elements() == false
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).

Stacktrace:
---
#0  ./MF_FES_RK4-Test: void dealii::internal::VectorDataExchange<2, double, 
dealii::VectorizedArray 
>::compress_start, 
(dealii::LinearAlgebra::distributed::Vector*)0>(unsigned int, 
dealii::LinearAlgebra::distributed::Vector&)
#1  ./MF_FES_RK4-Test: void dealii::internal::compress_start<2, 
dealii::LinearAlgebra::distributed::Vector, double, dealii::VectorizedArray, 
(dealii::LinearAlgebra::distributed::Vector*)0>(dealii::LinearAlgebra::distributed::Vector&, dealii::internal::VectorDataExchange<2, 
double, dealii::VectorizedArray >&, unsigned int)
#2  ./MF_FES_RK4-Test: dealii::internal::MFWorker >, 
dealii::LinearAlgebra::distributed::Vector, 
dealii::LinearAlgebra::distributed::Vector, Step40::LaplaceOperator<2, 2, 4>, 
true>::vector_compress_start()
#3  /opt/dealii/lib/libdeal_II.g.so.9.2.0-pre: 
dealii::internal::MatrixFreeFunctions::MPICommunication::execute()
#4 
 
/opt/intel/compilers_and_libraries_2019.5.281/linux/tbb/lib/intel64_lin/gcc4.7/libtbb_debug.so.2:
 

#5 
 
/opt/intel/compilers_and_libraries_2019.5.281/linux/tbb/lib/intel64_lin/gcc4.7/libtbb_debug.so.2:
 

#6 
 
/opt/intel/compilers_and_libraries_2019.5.281/linux/tbb/lib/intel64_lin/gcc4.7/libtbb_debug.so.2:
 

#7 
 
/opt/intel/compilers_and_libraries_2019.5.281/linux/tbb/lib/intel64_lin/gcc4.7/libtbb_debug.so.2:
 

#8 
 
/opt/intel/compilers_and_libraries_2019.5.281/linux/tbb/lib/intel64_lin/gcc4.7/libtbb_debug.so.2:
 

#9  /lib64/libpthread.so.0: 
#10  /lib64/libc.so.6: clone

Why that, and how can I fix it? 
Thanks!

-- 
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/48a874c7-af99-4634-bcf0-68e385be7cd6%40googlegroups.com.

[deal.II] Implementation of non-linear diffusion equation using the TimeStepping-Class and matrix-free routines

['Maxi Miller' via deal.II User Group]

I tried to implement a solver for the non-linear diffusion equation 
(\partial_t u = grad(u(grad u)) - f) using the TimeStepping-Class, the 
EmbeddedExplicitRungeKutta-Method and (for assembly) the matrix-free 
approach. For initial tests I used the linear heat equation with the 
solution u = exp(-pi * pi * t) * sin(pi * x) * sin(pi * y) and the assembly 
routine
template 
void LaplaceOperator::local_apply_cell(
const MatrixFree &   data,
vector_t &  dst,
const vector_t ,
const std::pair & cell_range) 
const
{
FEEvaluation 
phi(data);

for (unsigned int cell = cell_range.first; cell < cell_range.second; 
++cell)
{
phi.reinit(cell);
phi.read_dof_values_plain(src);
phi.evaluate(false, true);
for (unsigned int q = 0; q < phi.n_q_points; ++q)
{
auto gradient_coefficient = calculate_gradient_coefficient(
phi.get_gradient(q));
phi.submit_gradient(gradient_coefficient, q);
}
phi.integrate_scatter(false, true, dst);
}
}


and 
template 
inline DEAL_II_ALWAYS_INLINE
Tensor<1, n_components_to_use, Tensor<1, dim, Number>> 
calculate_gradient_coefficient(
#if defined(USE_NONLINEAR) || defined(USE_ADVECTION)
const Tensor<1, n_components_to_use, Number> _value,
#endif
const Tensor<1, n_components_to_use, Tensor<1, dim, Number>> &
input_gradient){
Tensor<1, n_components_to_use, Tensor<1, dim, Number>> ret_val;
for(size_t component = 0; component < n_components_to_use; ++
component){
for(size_t d = 0; d < dim; ++d){
ret_val[component][d] = -1. / (2 * M_PI * M_PI) * 
input_gradient[component][d];
}
}
return ret_val;
}


This approach works, and delivers correct results. Now I wanted to test the 
same approach for the non-linear diffusion equation with f = -exp(-2 * pi^2 
* t) * 0.5 * pi^2 * (-cos(2 * pi * (x - y)) - cos(2 * pi * (x + y)) + cos(2 
* pi * x) + cos(2 * pi * y)), which should be the solution to grad(u (grad 
u)) with u = exp(-pi^2*t) * sin(pi * x) * sin(pi * y). Thus, I changed the 
routines to
template 
void LaplaceOperator::local_apply_cell(
const MatrixFree &   data,
vector_t &  dst,
const vector_t ,
const std::pair & cell_range) 
const
{
FEEvaluation 
phi(data);

for (unsigned int cell = cell_range.first; cell < cell_range.second; 
++cell)
{
phi.reinit(cell);
phi.read_dof_values_plain(src);
phi.evaluate(true, true);
for(size_t q = 0; q < phi.n_q_points; ++q){
auto value = phi.get_value(q);
auto gradient = phi.get_gradient(q);
phi.submit_value(calculate_value_coefficient(value,
 phi.
quadrature_point(q),
 local_time), q
);
phi.submit_gradient(calculate_gradient_coefficient(value,
   gradient
), q);
}
phi.integrate_scatter(true, true, dst);
}
}

and 
template 
inline DEAL_II_ALWAYS_INLINE
Tensor<1, n_components_to_use, Number> calculate_value_coefficient(const 
Tensor<1, n_components_to_use, Number> _value,
   const 
Point ,
   const 
double ){
Tensor<1, n_components_to_use, Number> ret_val;
//(void) input_value;
(void) input_value;
for(size_t component = 0; component < n_components_to_use; ++
component){
const double x = point[component][0];
const double y = point[component][1];
ret_val[component] = (- exp(-2 * M_PI * M_PI * time)
  * 0.5 * M_PI * M_PI * (-cos(2 * M_PI * (x 
- y))
 - cos(2 * M_PI * (x 
+ y))
 + cos(2 * M_PI * x)
 + cos(2 * M_PI * y
)))
;
}
return ret_val;
}

template 
inline DEAL_II_ALWAYS_INLINE
Tensor<1, n_components_to_use, Tensor<1, dim, Number>> 
calculate_gradient_coefficient(
#if defined(USE_NONLINEAR) || defined(USE_ADVECTION)
const Tensor<1, n_components_to_use, Number> _value,
#endif
const Tensor<1, n_components_to_use, Tensor<1, dim, Number>> &
input_gradient){
Tensor<1, n_components_to_use, Tensor<1, dim, Number>> ret_val;
for(size_t component = 0; component < 

[deal.II] Implementation of a solver for the non-linear diffusion equation using the TimeStepping-Class

['Maxi Miller' via deal.II User Group]

I tried to implement a solver for the non-linear diffusion equation 
(\partial_t u = grad(u(grad u)) - f) using the TimeStepping-Class, the 
EmbeddedExplicitRungeKutta-Method and (for assembly) the matrix-free 
approach. For initial tests I used the linear heat equation with the 
solution u = exp(-pi * pi * t) * sin(pi * x) * sin(pi * y) and the assembly 
routine
template 
void LaplaceOperator::local_apply_cell(
const MatrixFree &   data,
vector_t &  dst,
const vector_t ,
const std::pair & cell_range) 
const
{
FEEvaluation 
phi(data);

for (unsigned int cell = cell_range.first; cell < cell_range.second; 
++cell)
{
phi.reinit(cell);
phi.read_dof_values_plain(src);
phi.evaluate(false, true);
for (unsigned int q = 0; q < phi.n_q_points; ++q)
{
auto gradient_coefficient = calculate_gradient_coefficient(
phi.get_gradient(q));
phi.submit_gradient(gradient_coefficient, q);
}
phi.integrate_scatter(false, true, dst);
}
}


and 
template 
inline DEAL_II_ALWAYS_INLINE
Tensor<1, n_components_to_use, Tensor<1, dim, Number>> 
calculate_gradient_coefficient(
#if defined(USE_NONLINEAR) || defined(USE_ADVECTION)
const Tensor<1, n_components_to_use, Number> _value,
#endif
const Tensor<1, n_components_to_use, Tensor<1, dim, Number>> &
input_gradient){
Tensor<1, n_components_to_use, Tensor<1, dim, Number>> ret_val;
for(size_t component = 0; component < n_components_to_use; ++
component){
for(size_t d = 0; d < dim; ++d){
ret_val[component][d] = -1. / (2 * M_PI * M_PI) * 
input_gradient[component][d];
}
}
return ret_val;
}


This approach works, and delivers correct results. Now I wanted to test the 
same approach for the non-linear diffusion equation with f = -exp(-2 * pi^2 
* t) * 0.5 * pi^2 * (-cos(2 * pi * (x - y)) - cos(2 * pi * (x + y)) + cos(2 
* pi * x) + cos(2 * pi * y)), which should be the solution to grad(u (grad 
u)) with u = exp(-pi^2*t) * sin(pi * x) * sin(pi * y). Thus, I changed the 
routines to
template 
void LaplaceOperator::local_apply_cell(
const MatrixFree &   data,
vector_t &  dst,
const vector_t ,
const std::pair & cell_range) 
const
{
FEEvaluation 
phi(data);

for (unsigned int cell = cell_range.first; cell < cell_range.second; 
++cell)
{
phi.reinit(cell);
phi.read_dof_values_plain(src);
phi.evaluate(true, true);
for(size_t q = 0; q < phi.n_q_points; ++q){
auto value = phi.get_value(q);
auto gradient = phi.get_gradient(q);
phi.submit_value(calculate_value_coefficient(value,
 phi.
quadrature_point(q),
 local_time), q
);
phi.submit_gradient(calculate_gradient_coefficient(value,
   gradient
), q);
}
phi.integrate_scatter(true, true, dst);
}
}

and 
template 
inline DEAL_II_ALWAYS_INLINE
Tensor<1, n_components_to_use, Number> calculate_value_coefficient(const 
Tensor<1, n_components_to_use, Number> _value,
   const 
Point ,
   const 
double ){
Tensor<1, n_components_to_use, Number> ret_val;
//(void) input_value;
(void) input_value;
for(size_t component = 0; component < n_components_to_use; ++
component){
const double x = point[component][0];
const double y = point[component][1];
ret_val[component] = (- exp(-2 * M_PI * M_PI * time)
  * 0.5 * M_PI * M_PI * (-cos(2 * M_PI * (x 
- y))
 - cos(2 * M_PI * (x 
+ y))
 + cos(2 * M_PI * x)
 + cos(2 * M_PI * y
)))
;
}
return ret_val;
}

template 
inline DEAL_II_ALWAYS_INLINE
Tensor<1, n_components_to_use, Tensor<1, dim, Number>> 
calculate_gradient_coefficient(
#if defined(USE_NONLINEAR) || defined(USE_ADVECTION)
const Tensor<1, n_components_to_use, Number> _value,
#endif
const Tensor<1, n_components_to_use, Tensor<1, dim, Number>> &
input_gradient){
Tensor<1, n_components_to_use, Tensor<1, dim, Number>> ret_val;
for(size_t component = 0; component < 

Re: [deal.II] Segfault when calling gather_evaluate() or read_dof_values_plain()

['Maxi Miller' via deal.II User Group]

For approximating the inverse mass matrix I followed step-48 with some 
modifications:
In the reinit-subroutine I added 
data.initialize_dof_vector(inv_mass_matrix);
FEEvaluation fe_eval(data);
const unsigned int   n_q_points = fe_eval.n_q_points;
for (unsigned int cell = 0; cell < data.n_cell_batches(); ++cell)
{
fe_eval.reinit(cell);
for (unsigned int q = 0; q < n_q_points; ++q)
fe_eval.submit_value(make_vectorized_array(1.), q);
fe_eval.integrate(true, false);
fe_eval.distribute_local_to_global(inv_mass_matrix);
}
inv_mass_matrix.compress(VectorOperation::add);
for (unsigned int k = 0; k < inv_mass_matrix.local_size(); ++k)
if (inv_mass_matrix.local_element(k) > 1e-15)
inv_mass_matrix.local_element(k) =
1. / inv_mass_matrix.local_element(k);
else
inv_mass_matrix.local_element(k) = 1;

and for applying the matrix I wrote
template 
void LaplaceOperator::local_apply_inverse_mass_matrix(
const MatrixFree &   data,
LinearAlgebra::distributed::Vector &  dst,
const LinearAlgebra::distributed::Vector ,
const std::pair & cell_range) 
const
{
dst.reinit(src);
dst = src;
dst.scale(inv_mass_matrix);
}
so that I could reuse the cell_loops from the original code. Still, the 
results indicate that something goes wrong. Should I handle the inverse 
matrix somehow different?
Thanks!
Am Dienstag, 14. Januar 2020 12:43:40 UTC+1 schrieb Martin Kronbichler:
>
> Dear Maxi,
>
> It looks like you are using a scalar finite element (FE_Q), but you 
> set FEEvaluation to operate on a vectorial field, i.e., 
>
> FEEvaluation phi(data);
>
> Notice the 4th template parameter "dim", which should be one. I agree it 
> is unfortunate that we do not provide a better error message, so I opened 
> an issue for it, https://github.com/dealii/dealii/issues/9312
>
> If you switch the parameter to 1, it should work. However, I want to point 
> out that you use a continuous element with CellwiseInverseMassMatrix - this 
> will not give you a valid matrix inverse. You need to either use a diagonal 
> mass matrix (see step-48) or use the cellwise inverse as a preconditioner 
> (when combined with suitable scaling) and solve the mass matrix iteratively.
>
> Best,
> Martin
> On 14.01.20 10:49, 'Maxi Miller' via deal.II User Group wrote:
>
> I could narrow it down to the function
>   
>  // same as above for has_partitioners_are_compatible == true
>   template <
> typename VectorType,
> typename std::enable_if::
> value,
> VectorType>::type * = nullptr>
>   inline void
>   check_vector_compatibility(
> const VectorType &vec,
> const internal::MatrixFreeFunctions::DoFInfo _info)
>   {
> (void)vec;
> (void)dof_info;
> Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
>ExcMessage(
>  "The parallel layout of the given vector is not "
>  "compatible with the parallel partitioning in MatrixFree. "
>  "Use MatrixFree::initialize_dof_vector to get a "
>  "compatible vector."));
>   }
>
> located in vector_access_internal.h. Apparently I have a segfault in the 
> function 
> 
>  template 
> inline bool
> Vector::partitioners_are_compatible(
>   const Utilities::MPI::Partitioner ) const
> {
>   return partitioner->is_compatible(part);
> }
>
> located in la_parallel_vector.templates.h, and thereby directly stopping 
> the program without triggering the assert()-function. The direct gdb output 
> for this function is 
> #3  0x700c0c5c in 
> dealii::LinearAlgebra::distributed::Vector dealii::MemorySpace::Host>::partitioners_are_compatible (this=0x0, 
> part=...) at 
> ~Downloads/git-files/dealii/include/deal.II/lac/la_parallel_vector.templates.h:2021
>
> I am not sure what that means. Could it point to a nullpointer somewhere 
> (after this=0x0)? Though, when putting a breakpoint to the first function 
> and printing the involved variables there, I get
> (gdb) print dof_info.vector_partitioner
> $2 = std::shared_ptr (use 
> count 3, weak count 0) = {get() = 0x7a60c0}
> (gdb) print vec.partitioner
> $3 = std::shared_ptr (use 
> count 3, weak count 0) = {get() = 0x7a60c0}
>
> i.e. no null ptr. Could the problem be there, or somewhere else?
> Thanks!
>
> Am Montag, 13. Januar 2020 21:41:54 UTC+1 schrieb Wolfgang Ban

Re: [deal.II] Segfault when calling gather_evaluate() or read_dof_values_plain()

['Maxi Miller' via deal.II User Group]

I could narrow it down to the function
 
 // same as above for has_partitioners_are_compatible == true
  template <
typename VectorType,
typename std::enable_if::
value,
VectorType>::type * = nullptr>
  inline void
  check_vector_compatibility(
const VectorType &vec,
const internal::MatrixFreeFunctions::DoFInfo _info)
  {
(void)vec;
(void)dof_info;
Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
   ExcMessage(
 "The parallel layout of the given vector is not "
 "compatible with the parallel partitioning in MatrixFree. "
 "Use MatrixFree::initialize_dof_vector to get a "
 "compatible vector."));
  }

located in vector_access_internal.h. Apparently I have a segfault in the 
function 
   
 template 
inline bool
Vector::partitioners_are_compatible(
  const Utilities::MPI::Partitioner ) const
{
  return partitioner->is_compatible(part);
}

located in la_parallel_vector.templates.h, and thereby directly stopping 
the program without triggering the assert()-function. The direct gdb output 
for this function is 
#3  0x700c0c5c in 
dealii::LinearAlgebra::distributed::Vector::partitioners_are_compatible (this=0x0, 
part=...) at 
~Downloads/git-files/dealii/include/deal.II/lac/la_parallel_vector.templates.h:2021

I am not sure what that means. Could it point to a nullpointer somewhere 
(after this=0x0)? Though, when putting a breakpoint to the first function 
and printing the involved variables there, I get
(gdb) print dof_info.vector_partitioner
$2 = std::shared_ptr (use count 3
, weak count 0) = {get() = 0x7a60c0}
(gdb) print vec.partitioner
$3 = std::shared_ptr (use count 3
, weak count 0) = {get() = 0x7a60c0}

i.e. no null ptr. Could the problem be there, or somewhere else?
Thanks!

Am Montag, 13. Januar 2020 21:41:54 UTC+1 schrieb Wolfgang Bangerth:
>
> On 1/13/20 8:41 AM, 'Maxi Miller' via deal.II User Group wrote: 
> > Therefore, I assume that I have some bugs in my code, but am not 
> experienced 
> > enough in writing matrix-free code for being able to debug it myself. 
> What 
> > kind of approach could I do now, or is there simply a glaring bug which 
> I did 
> > not see yet? 
>
> I haven't even looked at the code yet, but for segfaults there is a 
> relatively 
> simple cure: Run your code in a debugger. A segmentation fault is a 
> problem 
> where some piece of code tries to access data at an address (typically 
> through 
> a pointer) that is not readable or writable. The general solution to 
> finding 
> out what the issue is is to run the code in a debugger -- the debugger 
> will 
> then stop at the place where the problem happens, and you can inspect all 
> of 
> the values of the pointers and variables that are live at that location. 
> Once 
> you see what variable is at fault, the problem is either already obvious, 
> or 
> you can start to think about how a pointer got the value it has and why. 
>
> 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/d55e53b2-a0c8-4cc0-ab84-7521f8b23d81%40googlegroups.com.

[deal.II] Segfault when calling gather_evaluate() or read_dof_values_plain()

['Maxi Miller' via deal.II User Group]

I tried to implement a RK4-solver similar to the upcoming example 67 using 
the matrix-free technique, and compare it to other solvers I wrote earlier. 
For this purpose I chose the classical heat-equation for testing. When 
running the attached example, though, I get a segfault when calling either 
gather_evaluate() or read_dof_values_plain(). By replacing the code from 
read_dof_values_plain() with the code from the include file, I could narrow 
the problem down to (presumably) the operation read_write_operation(). When 
commenting this function out, the code continues until evaluate(), when 
not, it segfaults.
Therefore, I assume that I have some bugs in my code, but am not 
experienced enough in writing matrix-free code for being able to debug it 
myself. What kind of approach could I do now, or is there simply a glaring 
bug which I did not see yet?
Thanks!

-- 
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/06afdccd-8520-435a-9c49-2cfe1021998e%40googlegroups.com.
/* -
 *
 * Copyright (C) 2009 - 2019 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, redistribute
 * it, and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * The full text of the license can be found in the file LICENSE.md at
 * the top level directory of deal.II.
 *
 * -
 *
 * Author: Wolfgang Bangerth, Texas A University, 2009, 2010
 * Timo Heister, University of Goettingen, 2009, 2010
 */
#include 
#include 
#include 
#include 

#include 

#include 
#include 
#include 

namespace LA
{
#undef DEAL_II_WITH_PETSC
#if defined(DEAL_II_WITH_PETSC) && !defined(DEAL_II_PETSC_WITH_COMPLEX) && \
	!(defined(DEAL_II_WITH_TRILINOS) && defined(FORCE_USE_OF_TRILINOS))
	using namespace dealii::LinearAlgebraPETSc;
#  define USE_PETSC_LA
#elif defined(DEAL_II_WITH_TRILINOS)
	using namespace dealii::LinearAlgebraTrilinos;
#else
#  error DEAL_II_WITH_PETSC or DEAL_II_WITH_TRILINOS required
#endif
} // namespace LA
#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 

#include 
#include 

#include 
#include 

constexpr double time_step = 0.01;

constexpr unsigned int fe_degree= 2;
constexpr unsigned int n_q_points_1d= fe_degree + 2;
using Number= double;

enum LowStorageRungeKuttaScheme
{
	stage_3_order_3, /* Kennedy, Carpenter, Lewis, 2000 */
	stage_5_order_4, /* Kennedy, Carpenter, Lewis, 2000 */
	stage_7_order_4, /* Tselios, Simos, 2007 */
	stage_9_order_5, /* Kennedy, Carpenter, Lewis, 2000 */
};
constexpr LowStorageRungeKuttaScheme lsrk_scheme = stage_5_order_4;

namespace Step40
{
	using namespace dealii;

	class LowStorageRungeKuttaIntegrator
	{
	public:
		LowStorageRungeKuttaIntegrator(const LowStorageRungeKuttaScheme scheme)
		{
			if (scheme == stage_3_order_3)
			{
bi = {{0.245170287303492, 0.184896052186740, 0.569933660509768}};
ai = {{0.755726351946097, 0.386954477304099}};
			}
			else if (scheme == stage_5_order_4)
			{
bi = {{1153189308089. / 22510343858157.,
	   1772645290293. / 4653164025191.,
	   -1672844663538. / 4480602732383.,
	   2114624349019. / 3568978502595.,
	   5198255086312. / 14908931495163.}};
ai = {{970286171893. / 4311952581923.,
	   6584761158862. / 12103376702013.,
	   2251764453980. / 15575788980749.,
	   26877169314380. / 34165994151039.}};
			}
			else if (scheme == stage_7_order_4)
			{
bi = {{0.0941840925477795334,
	   0.149683694803496998,
	   0.285204742060440058,
	   -0.122201846148053668,
	   0.0605151571191401122,
	   0.345986987898399296,
	   0.186627171718797670}};
ai = {{0.241566650129646868 + bi[0],
	   0.0423866513027719953 + bi[1],
	   0.215602732678803776 + bi[2],
	   0.232328007537583987 + bi[3],
	   0.256223412574146438 + bi[4],
	   0.0978694102142697230 + bi[5]}};
			}
			else if (scheme == stage_9_order_5)
			{
bi = {{2274579626619. / 23610510767302.,
	   693987741272. / 12394497460941.,
	   -347131529483. / 15096185902911.,
	   1144057200723. / 32081666971178.,
	   

[deal.II] Re: "fatal error: MueLu_CreateEpetraPreconditioner.hpp: No such file or directory" when compiling with Trilinos

['Maxi Miller' via deal.II User Group]

I ran into several other compilation problems with that trilinos version, 
thus I decided to try the suggested patch instead.

Am Donnerstag, 9. Januar 2020 16:31:21 UTC+1 schrieb Bruno Turcksin:
>
> Hi,
>
> This is a known problem with the development version of Trilinos 
> https://github.com/dealii/dealii/pull/9237 Use the latest release of 
> Trilinos instead.
>
> Best,
>
> Bruno
>
> On Thursday, January 9, 2020 at 10:06:13 AM UTC-5, Maxi Miller wrote:
>>
>> Hei,
>> when trying to compile deal.II with one of the latest Trilinos-Versions, 
>> compilation fails with "fatal error: MueLu_CreateEpetraPreconditioner.hpp: 
>> No such file or directory", but the configuration process goes through 
>> without rising any problems. Which modules do I have to enable in Trilinos 
>> to fix this warning?
>> Thanks!
>>
>

-- 
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/cb7e1a0d-819c-484b-9097-404fbcef2b44%40googlegroups.com.

[deal.II] "fatal error: MueLu_CreateEpetraPreconditioner.hpp: No such file or directory" when compiling with Trilinos

['Maxi Miller' via deal.II User Group]

Hei,
when trying to compile deal.II with one of the latest Trilinos-Versions, 
compilation fails with "fatal error: MueLu_CreateEpetraPreconditioner.hpp: 
No such file or directory", but the configuration process goes through 
without rising any problems. Which modules do I have to enable in Trilinos 
to fix this warning?
Thanks!

-- 
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/f904a911-77e1-48cb-a60f-fc95499e7668%40googlegroups.com.

[deal.II] Usage of the AD framework for solving a time-dependent equation

['Maxi Miller' via deal.II User Group]

I tried to use the AD framework in my own codes (mainly for solving 
nonlinear equations), and therefore wrote a small test program for solving 
the linear diffusion equation (after I can calculate the jacobian by hand 
for that equation). Thus, my residual for the AD-jacobian is -0.5 * (\nabla 
u_1 + \nabla u_0)\cdot\nabla\varphi - \frac{1}{dt}(u_1 - u_0)\cdot\varphi, 
with u_0 my current solution, and u_1 my solution for the next step. When 
creating the jacobian by hand, I use the same approach as shown in example 
33, and assume that u_1 = u_0 in my first step, resulting in 
-0.5\cdot(\nabla\varphi_i\cdot\nabla\varphi_j) - 
\frac{1}{dt}\cdot\varphi_i\cdot\varphi_j, and my right hand side is \nabla 
u_0\cdot\nabla\varphi.
Solving both systems should give me a value for \delta u, i.e. the update 
for my current solution, such that u_1 = \delta u + u_0 (after the system 
is linear).
Though, when running the code in the attachment, I get the correct result 
(i.e. \delta u) for my hand-assembled system, but the solution I get for 
the AD-generated jacobian is u_1, not \delta u. Did I misinterpret 
something, or is there another bug in the code?
Thanks!

-- 
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/e556dfdd-d5e8-453c-896a-0fb2f3cc73b2%40googlegroups.com.
/* -
 *
 * Copyright (C) 2009 - 2019 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, redistribute
 * it, and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * The full text of the license can be found in the file LICENSE.md at
 * the top level directory of deal.II.
 *
 * -
 *
 * Author: Wolfgang Bangerth, Texas A University, 2009, 2010
 * Timo Heister, University of Goettingen, 2009, 2010
 */
#include 
#include 
#include 
#include 

namespace LA
{
#undef DEAL_II_WITH_PETSC
#if defined(DEAL_II_WITH_PETSC) && !defined(DEAL_II_PETSC_WITH_COMPLEX) && \
	!(defined(DEAL_II_WITH_TRILINOS) && defined(FORCE_USE_OF_TRILINOS))
	using namespace dealii::LinearAlgebraPETSc;
#  define USE_PETSC_LA
#elif defined(DEAL_II_WITH_TRILINOS)
	using namespace dealii::LinearAlgebraTrilinos;
#else
#  error DEAL_II_WITH_PETSC or DEAL_II_WITH_TRILINOS required
#endif
} // namespace LA
#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 
#include 

#include 
#include 

constexpr double time_step = 0.0001;

namespace Step40
{
	using namespace dealii;


	template 
	class InitialCondition : public Function
	{
	public:
		InitialCondition() : Function(1){

		}
		virtual double value(const Point , const unsigned int component) const override;
		virtual Tensor<1, dim> gradient(const Point , const unsigned int component) const override;
	};

	template 
	double InitialCondition::value(const Point , const unsigned int) const{
		const double x = p[0];
		const double y = p[1];

		return sin(M_PI * x) * sin(M_PI * y);
	}

	template 
	Tensor<1, dim> InitialCondition::gradient(const Point , const unsigned int) const{
		const double x = p[0];
		const double y = p[1];

		Tensor<1, dim> return_value;
		return_value[0] = M_PI * cos(M_PI * x) * sin(M_PI * y);
		return_value[1] = M_PI * sin(M_PI * x) * cos(M_PI * y);

		return return_value;
	}

	template 
	class LaplaceProblem
	{
	public:
		LaplaceProblem();
		void run();
	private:
		void setup_system();
		void assemble_jacobian_by_hand_dynamic();
		void assemble_jacobian_by_AD_dynamic();
		void solve();
		void solve_with_newton_hand();
		void solve_with_newton_AD();
		void refine_grid();
		void output_results(const unsigned int cycle) const;

		MPI_Comm mpi_communicator;
		parallel::distributed::Triangulation triangulation;
		FE_Q   fe;
		DoFHandler dof_handler;
		IndexSet locally_owned_dofs;
		IndexSet locally_relevant_dofs;
		AffineConstraints constraints;
		LA::MPI::SparseMatrix system_matrix, jacobian_matrix_AD, jacobian_matrix_hand, mass_matrix;
		LA::MPI::Vector   locally_relevant_solution,
		locally_relevant_solution_AD,
		locally_relevant_solution_hand;
		LA::MPI::Vector   system_rhs, residual_rhs_AD, residual_rhs_hand;
		ConditionalOStream pcout;
		TimerOutput   

Re: [deal.II] Arpack solver reports 4294967295 iterations before converging

['Maxi Miller' via deal.II User Group]

Like this: https://github.com/dealii/dealii/pull/8899?

Am Freitag, 11. Oktober 2019 16:42:24 UTC+2 schrieb Wolfgang Bangerth:
>
> On 10/11/19 7:22 AM, 'Maxi Miller' via deal.II User Group wrote: 
> > It is using the SolverControl-class for setting parameters such as the 
> > maximum iterations, but the used iterations are never returned (no call 
> > to check()). Therefore, was there a reason for setting the 
> > control()-object as constant, or can I change it to a non-constant 
> value? 
>
> I think we probably didn't think this through. Do you happen to know 
> whether ARPACK returns the number of iterations and final 
> residual/error? If so, then what we should probably call the check() 
> function with these values, ignoring the return code, just so that the 
> values are stored in the SolverControl object. 
>
> Want to give that a try and make it into a patch? 
>
> 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/173b293b-76fe-4f5f-aa4c-9357172a4429%40googlegroups.com.

Re: [deal.II] Arpack solver reports 4294967295 iterations before converging

['Maxi Miller' via deal.II User Group]

It is using the SolverControl-class for setting parameters such as the 
maximum iterations, but the used iterations are never returned (no call to 
check()). Therefore, was there a reason for setting the control()-object as 
constant, or can I change it to a non-constant value?

Am Freitag, 11. Oktober 2019 07:35:41 UTC+2 schrieb Wolfgang Bangerth:
>
> On 10/9/19 7:47 AM, 'Maxi Miller' via deal.II User Group wrote: 
> > This number already comes from solver_control, given as 
> > 
> > solver_control.last_step (); 
>
> Hm, that's really quite unfortunate. It took me a couple minutes of 
> looking 
> into the code, but I *think* what is happening is that the ArpackSolver 
> class 
> simply doesn't use the solver control object at all, and that the -1 for 
> the 
> number of iterations is simply the (unchanged) initializer the 
> SolverControl 
> class uses for the number of iterations. So that would mean that it 
> doesn't 
> actually indicate an error. 
>
> Can you check whether my interpretation is correct? If so, we should 
> probably 
> just remove the argument from the ArpackSolver class's constructor, as 
> well as 
> the corresponding member variable. 
>
> 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/6221cc43-db93-47dc-b2de-6e237b03d6e4%40googlegroups.com.

Re: [deal.II] Arpack solver reports 4294967295 iterations before converging

['Maxi Miller' via deal.II User Group]

This number already comes from solver_control, given as 

solver_control.last_step ();


Am Montag, 7. Oktober 2019 15:33:29 UTC+2 schrieb Wolfgang Bangerth:
>
> On 10/7/19 5:32 AM, 'Maxi Miller' via deal.II User Group wrote: 
> > I implemented the changes in step-36 suggested in the tutorial for using 
> > ARPACK instead of PETSc (as shown in the attachment). Now, the solver 
> reports 
> > a convergence in 4294967295 iterations (with correct result), thus I 
> assume a 
> > overflow/underflow bug. Is that something which I should report in the 
> github 
> > issues? 
>
> 4294967295=-1 when expressed in unsigned int. This usually indicates some 
> kind 
> of error condition. It would be worth figuring out where this number is 
> generated because it's pretty clear that that is not the actual number of 
> iterations performed. 
>
> 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/3ba1cfa0-7a23-482a-8a47-ae24668d2df1%40googlegroups.com.

[deal.II] Arpack solver reports 4294967295 iterations before converging

['Maxi Miller' via deal.II User Group]

I implemented the changes in step-36 suggested in the tutorial for using 
ARPACK instead of PETSc (as shown in the attachment). Now, the solver 
reports a convergence in 4294967295 iterations (with correct result), thus 
I assume a overflow/underflow bug. Is that something which I should report 
in the github issues?
Thanks!

-- 
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/c71acebf-2574-4bff-bb24-e34be30bcaa4%40googlegroups.com.
/* -
 *
 * Copyright (C) 2009 - 2019 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, redistribute
 * it, and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * The full text of the license can be found in the file LICENSE.md at
 * the top level directory of deal.II.
 *
 * -

 *
 * Authors: Toby D. Young, Polish Academy of Sciences,
 *  Wolfgang Bangerth, Texas A University
 */

// @sect3{Include files}

// As mentioned in the introduction, this program is essentially only a
// slightly revised version of step-4. As a consequence, most of the following
// include files are as used there, or at least as used already in previous
// tutorial programs:
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 

// IndexSet is used to set the size of each PETScWrappers::MPI::Vector:
#include 

// PETSc appears here because SLEPc depends on this library:
#include 
#include 

// And then we need to actually import the interfaces for solvers that SLEPc
// provides:
#include 

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

// We also need some standard C++:
#include 
#include 

// Finally, as in previous programs, we import all the deal.II class and
// function names into the namespace into which everything in this program
// will go:
namespace Step36
{
	using namespace dealii;

	// @sect3{The EigenvalueProblem class template}

	// Following is the class declaration for the main class template. It looks
	// pretty much exactly like what has already been shown in step-4:
	template 
	class EigenvalueProblem
	{
	public:
		EigenvalueProblem(const std::string _file);
		void run();

	private:
		void make_grid_and_dofs();
		void assemble_system();
		unsigned int solve();
		void output_results() const;

		Triangulation triangulation;
		FE_Q  fe;
		DoFHandlerdof_handler;

		// With these exceptions: For our eigenvalue problem, we need both a
		// stiffness matrix for the left hand side as well as a mass matrix for
		// the right hand side. We also need not just one solution function, but a
		// whole set of these for the eigenfunctions we want to compute, along
		// with the corresponding eigenvalues:
		//	PETScWrappers::SparseMatrix stiffness_matrix, mass_matrix;
		//	std::vector eigenfunctions;
		//	std::vector eigenvalues;
		SparsityPattern sparsity_pattern;
		SparseMatrixstiffness_matrix, mass_matrix;
		std::vector >eigenfunctions;
		std::vector>   eigenvalues;

		// And then we need an object that will store several run-time parameters
		// that we will specify in an input file:
		ParameterHandler parameters;

		// Finally, we will have an object that contains "constraints" on our
		// degrees of freedom. This could include hanging node constraints if we
		// had adaptively refined meshes (which we don't have in the current
		// program). Here, we will store the constraints for boundary nodes
		// $U_i=0$.
		AffineConstraints constraints;
	};

	// @sect3{Implementation of the EigenvalueProblem class}

	// @sect4{EigenvalueProblem::EigenvalueProblem}

	// First up, the constructor. The main new part is handling the run-time
	// input parameters. We need to declare their existence first, and then read
	// their values from the input file whose name is specified as an argument
	// to this function:
	template 
	EigenvalueProblem::EigenvalueProblem(const std::string _file)
		: fe(1)
		, dof_handler(triangulation)
	{
		// TODO investigate why the minimum number of refinement steps required to
		// obtain the correct eigenvalue degeneracies is 6
		

Re: [deal.II] Create a LinearAlgebra::distributed::Vector from a LinearAlgebra::distributed::Vector

['Maxi Miller' via deal.II User Group]

Yes, that worked, thanks!

Am Montag, 23. September 2019 17:17:09 UTC+2 schrieb Martin Kronbichler:
>
> Dear Maxi,
>
> Did you try to call 'operator=' in the same line as you declare the 
> variable? In that case, you still go through the (missing) constructor. It 
> should work if you first declare the float vector in one statement and then 
> assign in the next. We have this function:
>
>
> https://dealii.org/current/doxygen/deal.II/classLinearAlgebra_1_1distributed_1_1Vector.html#ad364c47c5b3b7b991490cd626c729f0b
>
> Furthermore, there is also a function 
> LinearAlgebra::distributed::Vector::copy_locally_owned_subrange_from() 
> which also works with different number types, say double to float:
>
>
> https://dealii.org/current/doxygen/deal.II/classLinearAlgebra_1_1distributed_1_1Vector.html#ab792ddb04b95a220e489f2d7f9eee990
>
> (In that case, no data is exchanged, and the vectors need to have 
> compatible locally owned dofs before calling.)
>
> Best,
> Martin
> On 23.09.19 17:07, 'Maxi Miller' via deal.II User Group wrote:
>
> When doing that in my code, I get the error
> error: conversion from ‘Vector’ to non-scalar type ‘Vector<
> float,[...]>’ requested
>
> thus I assume I have to implement it myself...
>
> Am Freitag, 13. September 2019 17:53:37 UTC+2 schrieb Wolfgang Bangerth: 
>>
>> On 9/13/19 8:39 AM, 'Maxi Miller' via deal.II User Group wrote: 
>> > In my program I have to reuse old solutions (due to time dependency), 
>> > but would like to use them both as double values and float values 
>> > (preconditioning does not require double values). Unfortunately my old 
>> > solution is only given as a 
>> > LinearAlgebra::distributed::Vector(), but I have to use a 
>> > LinearAlgebra::distributed::Vector. Direct conversion via 
>> > LinearAlgebra::distributed::Vector float_solution(old_solution) 
>> > does not work. Are there other ways for conversion? 
>>
>> I think that the other classes all have conversion constructors or 
>> conversion operator=. If this class doesn't have one, you probably want 
>> to add it to support what you need. It should be a relatively easy patch 
>> to write, and we'd of course be happy to take it, as always! 
>>
>> 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 dea...@googlegroups.com .
> To view this discussion on the web visit 
> https://groups.google.com/d/msgid/dealii/1e88472a-6dce-4a26-ada0-197a82adea63%40googlegroups.com
>  
> <https://groups.google.com/d/msgid/dealii/1e88472a-6dce-4a26-ada0-197a82adea63%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/07d2ac8c-96f4-4621-aca6-e29403f67a9b%40googlegroups.com.

Re: [deal.II] Create a LinearAlgebra::distributed::Vector from a LinearAlgebra::distributed::Vector

['Maxi Miller' via deal.II User Group]

When doing that in my code, I get the error
error: conversion from ‘Vector’ to non-scalar type ‘Vector<
float,[...]>’ requested

thus I assume I have to implement it myself...

Am Freitag, 13. September 2019 17:53:37 UTC+2 schrieb Wolfgang Bangerth:
>
> On 9/13/19 8:39 AM, 'Maxi Miller' via deal.II User Group wrote: 
> > In my program I have to reuse old solutions (due to time dependency), 
> > but would like to use them both as double values and float values 
> > (preconditioning does not require double values). Unfortunately my old 
> > solution is only given as a 
> > LinearAlgebra::distributed::Vector(), but I have to use a 
> > LinearAlgebra::distributed::Vector. Direct conversion via 
> > LinearAlgebra::distributed::Vector float_solution(old_solution) 
> > does not work. Are there other ways for conversion? 
>
> I think that the other classes all have conversion constructors or 
> conversion operator=. If this class doesn't have one, you probably want 
> to add it to support what you need. It should be a relatively easy patch 
> to write, and we'd of course be happy to take it, as always! 
>
> 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/1e88472a-6dce-4a26-ada0-197a82adea63%40googlegroups.com.

[deal.II] Create a LinearAlgebra::distributed::Vector from a LinearAlgebra::distributed::Vector

['Maxi Miller' via deal.II User Group]

In my program I have to reuse old solutions (due to time dependency), but 
would like to use them both as double values and float values 
(preconditioning does not require double values). Unfortunately my old 
solution is only given as a LinearAlgebra::distributed::Vector(), 
but I have to use a LinearAlgebra::distributed::Vector. Direct 
conversion via LinearAlgebra::distributed::Vector 
float_solution(old_solution) does not work. Are there other ways for 
conversion?

-- 
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/59bcb18e-354e-4747-a24a-424b67df3a2a%40googlegroups.com.

Re: [deal.II] Requirements for citing deal.II

['Maxi Miller' via deal.II User Group]

Ok, will do that, thanks!

Am Mittwoch, 11. September 2019 16:52:51 UTC+2 schrieb Wolfgang Bangerth:
>
> On 9/11/19 3:27 AM, 'Maxi Miller' via deal.II User Group wrote: 
> > In one of my side projects I include the library deal.II, but use it 
> > only for reading in parameter files and printing out data in the HDF5 
> > format. I included it mostly because it is in use in most of my main 
> > projects, and thus it already was available. Furthermore, being able to 
> > reuse parameter files and the rather simple way of storing data were the 
> > main reasons for including the library in this project. 
> > Thus, I was now wondering if that usage is enough to support citing the 
> > library in publications, even though I never use the main features of 
> it? 
>
> My general rule of thumb when I'm faced with similar dilemmas is "if I 
> use, I cite it". That's because someone still put it many weeks of work 
> to write the parameter handling and HDF5 capabilities (in your case; or 
> whatever else it is that I use), and it doesn't cost me very much to 
> cite some other paper. Indeed, my general approach (in this and any 
> other regard) is to be generous with praise. 
>
> 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/d92ed4cc-37dd-41c7-b498-28f16a9a9d4b%40googlegroups.com.

[deal.II] Requirements for citing deal.II

['Maxi Miller' via deal.II User Group]

In one of my side projects I include the library deal.II, but use it only 
for reading in parameter files and printing out data in the HDF5 format. I 
included it mostly because it is in use in most of my main projects, and 
thus it already was available. Furthermore, being able to reuse parameter 
files and the rather simple way of storing data were the main reasons for 
including the library in this project.
Thus, I was now wondering if that usage is enough to support citing the 
library in publications, even though I never use the main features of it?
Thanks!

-- 
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/051ece36-9853-403d-9c34-01e974413347%40googlegroups.com.

Re: [deal.II] VectorTools::Compress() fails after mesh_loop and distribute_local_to_global

['Maxi Miller' via deal.II User Group]

Will test that, just need to recompile my debug-version of deal.II

Am Samstag, 31. August 2019 18:17:34 UTC+2 schrieb Timo Heister:
>
> This is an interesting problem that might be a deal.II bug. I can't 
> see how it has anything to do with mesh_loop, but it looks like 
> operator=(0) is not setting last_action. 
>
> That said, I am confused because .reinit() should set last_action to 
> "Zero". 
>
> Can you go into the line of the assert in trilinos_vector.cc and print 
> the value of "last_action" on all processors? 
>
> If you can not figure it out, can you make a minimal complete example I 
> can try? 
>
>
> On Tue, Aug 27, 2019 at 6:59 AM 'Maxi Miller' via deal.II User Group 
> > wrote: 
> > 
> > I am trying to assemble my residual/RHS by using mesh_loop: 
> > evaluation_point = present_solution; 
> > extension_vector = newton_update; 
> > extension_vector *= alpha; 
> > 
> > evaluation_point += extension_vector; 
> > 
> > residual.reinit(locally_owned_dofs, mpi_communicator); 
> > residual_result.reinit(locally_owned_dofs, mpi_communicator); 
> > residual = 0; 
> > residual_result = 0; 
> > 
> > const QGauss  quadrature_formula(fe.degree + gauss_size); 
> > 
> > 
> > using CellFilter = FilteredIterator DoFHandler<2>::active_cell_iterator>; 
> > std::cout << "Looping over all cells\n"; 
> > 
> > MeshWorker::mesh_loop(dof_handler.begin_active(), 
> >   dof_handler.end(), 
> >   std::bind( ad_type_code>::local_assemble_residual, 
> > this, 
> > std::placeholders::_1, 
> > std::placeholders::_2, 
> > std::placeholders::_3), 
> >   std::bind( ad_type_code>::copy_local_to_global_residual, 
> > this, 
> > std::placeholders::_1), 
> >   ResidualScratchData(fe, 
> quadrature_formula, update_values | update_gradients | 
> update_quadrature_points | update_JxW_values, I_val, I_val_old), 
> >   ResidualCopyData()); 
> > 
> > std::cout << "Compressing residual\n"; 
> > residual.compress(VectorOperation::add); 
> > std::cout << "Residual compressed\n"; 
> > boundary_constraints.set_zero(residual); 
> > std::cout << "Boundary constraints set to zero\n"; 
> > residual_result = residual; 
> > 
> > 
> > 
> > with 
> > copy_local_to_global_residual 
> > 
> > defined as 
> > void MinimalSurfaceProblem ad_type_code>::copy_local_to_global_residual(const ResidualCopyData ) 
> > { 
> > 
> boundary_constraints.distribute_local_to_global(data.local_residual, 
> > 
> data.local_dof_indices, 
> > residual); 
> > } 
> > 
> > Unfortunately, my program fails after some iterations with the error 
> >  
> > An error occurred in line <624> of file 
> <~/Downloads/git-files/dealii/source/lac/trilinos_vector.cc> in function 
> > void 
> dealii::TrilinosWrappers::MPI::Vector::compress(dealii::VectorOperation::values)
>  
>
> > The violated condition was: 
> > result.max - result.min < 1e-5 
> > Additional information: 
> > Not all processors agree whether the last operation on this vector 
> was an addition or a set operation. This will prevent the compress() 
> operation from succeeding. 
> > 
> > Stacktrace: 
> > --- 
> > #0  /opt/dealii/lib/libdeal_II.g.so.9.2.0-pre: 
> dealii::TrilinosWrappers::MPI::Vector::compress(dealii::VectorOperation::values)
>  
>
> > #1  bin/TTM-equation: TTM_Calculation::MinimalSurfaceProblem<2, 
> Silicon_Verburg::physics_equations, 
> (dealii::Differentiation::AD::NumberTypes)3>::compute_residual(double, 
> dealii::TrilinosWrappers::MPI::Vector&) 
> > #2  bin/TTM-equation: TTM_Calculation::MinimalSurfaceProblem<2, 
> Silicon_Verburg::physics_equations, 
> (dealii::Differentiation::AD::NumberTypes)3>::recalculate_step_length() 
> > #3  bin/TTM-

[deal.II] VectorTools::Compress() fails after mesh_loop and distribute_local_to_global

['Maxi Miller' via deal.II User Group]

I am trying to assemble my residual/RHS by using mesh_loop:
evaluation_point = present_solution;
extension_vector = newton_update;
extension_vector *= alpha;

evaluation_point += extension_vector;

residual.reinit(locally_owned_dofs, mpi_communicator);
residual_result.reinit(locally_owned_dofs, mpi_communicator);
residual = 0;
residual_result = 0;

const QGauss  quadrature_formula(fe.degree + gauss_size);


using CellFilter = FilteredIterator::
active_cell_iterator>;
std::cout << "Looping over all cells\n";

MeshWorker::mesh_loop(dof_handler.begin_active(),
  dof_handler.end(),
  std::bind(::local_assemble_residual,
this,
std::placeholders::_1,
std::placeholders::_2,
std::placeholders::_3),
  std::bind(::copy_local_to_global_residual,
this,
std::placeholders::_1),
  ResidualScratchData(fe, quadrature_formula, 
update_values | update_gradients | update_quadrature_points | 
update_JxW_values, I_val, I_val_old),
  ResidualCopyData());

std::cout << "Compressing residual\n";
residual.compress(VectorOperation::add);
std::cout << "Residual compressed\n";
boundary_constraints.set_zero(residual);
std::cout << "Boundary constraints set to zero\n";
residual_result = residual;



with 
copy_local_to_global_residual

defined as 
void MinimalSurfaceProblem::
copy_local_to_global_residual(const ResidualCopyData )
{
boundary_constraints.distribute_local_to_global(data.local_residual,
data.
local_dof_indices,
residual);
}

Unfortunately, my program fails after some iterations with the error

An error occurred in line <624> of file <~/Downloads/git-files/dealii/source
/lac/trilinos_vector.cc> in function
void dealii::TrilinosWrappers::MPI::Vector::compress(dealii::
VectorOperation::values)
The violated condition was: 
result.max - result.min < 1e-5
Additional information: 
Not all processors agree whether the last operation on this vector was 
an addition or a set operation. This will prevent the compress() operation 
from succeeding.

Stacktrace:
---
#0  /opt/dealii/lib/libdeal_II.g.so.9.2.0-pre: 
dealii::TrilinosWrappers::MPI::Vector::compress(dealii::VectorOperation::values)
#1  bin/TTM-equation: TTM_Calculation::MinimalSurfaceProblem<2, 
Silicon_Verburg::physics_equations, 
(dealii::Differentiation::AD::NumberTypes)3>::compute_residual(double, 
dealii::TrilinosWrappers::MPI::Vector&)
#2  bin/TTM-equation: TTM_Calculation::MinimalSurfaceProblem<2, 
Silicon_Verburg::physics_equations, 
(dealii::Differentiation::AD::NumberTypes)3>::recalculate_step_length()
#3  bin/TTM-equation: TTM_Calculation::MinimalSurfaceProblem<2, 
Silicon_Verburg::physics_equations, 
(dealii::Differentiation::AD::NumberTypes)3>::update_newton_solution(double 
const&, double const&)
#4  bin/TTM-equation: TTM_Calculation::MinimalSurfaceProblem<2, 
Silicon_Verburg::physics_equations, 
(dealii::Differentiation::AD::NumberTypes)3>::run()
#5  bin/TTM-equation: main


Based on the output I assume that the distribute_local_to_global-function 
is the problem (when running with 8 threads, I get 8 lines with "Looping 
over all cells", but only 7 with "Compressing residual". The last operation 
(after doing a reinit() and setting it to 0) I apply onto the vector 
residual is distribute_local_to_global(), thus this is the current suspect. 
Or are there other other possibilities?

Note: I had a similar error earlier, but that time I forgot the final 
compress()-call before calling set_zero().

-- 
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/8301da2b-48e4-4600-8d37-a7a988d38a98%40googlegroups.com.

[deal.II] Mixed-precision algorithms in deal.II

['Maxi Miller' via deal.II User Group]

I tried to implement a simple mixed-precision algorithm in deal.II, which 
creates my matrix using float values, and my vectors using double values. 
Unfortunately, the function distribute_local_to_global(cell_matrix, 
cell_rhs, local_dof_indices, system_matrix, system_rhs) does not support 
mixed types. Is there an alternative function, or do I have to create two 
AffineConstraints-objects, initialize both with different values (double 
and float) and use them separately?

-- 
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/39fb6604-f0fc-43a7-b949-7e728e5ac14f%40googlegroups.com.

[deal.II] Creation of diagonal matrix for the heat equation in a matrix-free context

['Maxi Miller' via deal.II User Group]

Following step-37 I tried to extend it for non-linear equations, using a 
jacobian approximation for the "matrix" (which can be written as Jv = (F(u 
+ ev) - F(u)) * 1/e, with e a small value (~1e-6), F() the residual of the 
function, u the current solution, v the krylov vector obtained from the 
GMRES-solver and J the jacobian matrix). Using this method I tried to 
calculate the diagonal, which then is used for the preconditioner. 
Based on the equation given above and the equation for the linear heat 
equation, I get

F(u) = (u - u_0)/dt - (nabla^2u+nabla^2u_0)/2
F(u+ev) = F(u) + ev/dt - e(nabla^2v/2)

Thus Jv can be written as

Jv = (v/dt - nabla^2v/2)

or in weak form

Jv = v\phi/dt + 0.5*\nabla v\nabla\phi

Using that I tried to write the function for creating the diagonal 
accordingly:

template 
void JacobianOperator::local_compute_diagonal(
const MatrixFree & data,
vector_t ,
const unsigned int &,
const std::pair _range) const
{
FEEvaluation phi(data), 
phi_old(data);
AlignedVector> diagonal(phi.dofs_per_cell);
for (unsigned int cell = cell_range.first; cell < cell_range.second; ++
cell)
{
phi_old.reinit(cell);
for (unsigned int i = 0; i < phi.dofs_per_cell; ++i)
{
for (unsigned int j = 0; j < phi.dofs_per_cell; ++j){
phi_old.submit_dof_value(VectorizedArray(), j);
}
phi_old.submit_dof_value(make_vectorized_array(1.), i);
phi_old.evaluate(true, true);
for (unsigned int q = 0; q < phi_old.n_q_points; ++q){
phi_old.submit_gradient(0.5
* (phi_old.get_gradient(q)),
q);
phi_old.submit_value((phi_old.get_value(q) / time_step),
 q);
}
phi_old.integrate(true, true);
diagonal[i] = phi_old.get_dof_value(i);
}
for (unsigned int i = 0; i < phi_old.dofs_per_cell; ++i)
phi_old.submit_dof_value(diagonal[i], i);
phi_old.distribute_local_to_global(dst);
}
}

Still, the convergence behaviour is significantly worse than using an 
Identity preconditioner. Thus I was wondering if this is related to the 
wrong preconditioner, or due to wrong code? I am still trying to understand 
how to create the preconditioner for that case most efficiently, thus there 
might be some errors here.
Thanks!

-- 
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/cd41c84f-9441-4ca9-b48a-d8b9f120456f%40googlegroups.com.

Re: [deal.II] Re: MGTransferMatrixFree<> interpolate_to_mg() fails with segfault, when compiled in release mode

['Maxi Miller' via deal.II User Group]

Found my mistake, I forgot to call 
mg_transfer.build(dof_handler);
before doing anything with it.

Am Donnerstag, 1. August 2019 09:02:15 UTC+2 schrieb Maxi Miller:
>
> Hei,
> the program crashes in release mode way before it crashes in debug mode. 
> The code shown is completely from the deal.II-library, with some added 
> std::cout<<-lines. The only code written by myself is in main.cc, but the 
> crash happens in the code from the library.
>
>
> Am Donnerstag, 1. August 2019 05:27:21 UTC+2 schrieb Wolfgang Bangerth:
>>
>>
>> Maxi, 
>>
>> since your program is crashing in debug mode even before the location 
>> where it 
>> crashes in release mode, let's just ignore the release mode issue -- if a 
>> program doesn't run in debug mode, you should not assume that it does 
>> anything 
>> reasonable at all in release mode. 
>>
>>
>> > When compiling in debug mode, my results vary, and the program behaves 
>> > unpredictable, even though I run it in serial mode. In this case it 
>> fails 
>> > earlier during 
>> > 
>> > | 
>> >  std::cout <<"Looped over all levels, performing plain_copy: 
>> > "<> > "<> > if(perform_plain_copy) 
>> > { 
>> > // In this case, we can simply copy the local range. 
>> > AssertDimension(dst[dst.max_level()].local_size(),src.local_size()); 
>> >  dst[dst.max_level()].copy_locally_owned_data_from(src); 
>> > return; 
>> > } 
>> > elseif(perform_renumbered_plain_copy) 
>> > { 
>> > //std::cout << "dst_level-size: " << 
>> dst[dst.max_level()].local_size() 
>> > << '\n'; 
>> > //std::cout << "src-size: " << src.local_size() << '\n'; 
>> > Assert(copy_indices_level_mine.back().empty(),ExcInternalError()); 
>> > LinearAlgebra::distributed::Vector_level = 
>> >  dst[dst.max_level()]; 
>> > for(std::pairi :this_copy_indices.back()) 
>> >  
>> dst_level.local_element(i.second)=src.local_element(i.first); 
>> >  std::cout <<"Renumbered copy done\n"; 
>> > return; 
>> > } 
>> > | 
>> > 
>> > plain_copy is (at some point) equal to 0, resulting in the program 
>> choosing 
>> > the second option (perform_renumbered_plain_copy), where it crashes 
>> with a 
>> > segfault while evaluating the loop 
>> > | 
>> > 
>> > 
>> > for(std::pairi :this_copy_indices.back()) 
>> >  
>> dst_level.local_element(i.second)=src.local_element(i.first); 
>> > | 
>> > 
>> > This behaviour is rather random, and changes at different runs even 
>> though the 
>> > program was not recompiled. 
>> > 
>> > Why is that happening? Did I forget to initialize something, or is that 
>> an 
>> > internal bug? 
>>
>> Is the code you show somewhere in the library, or in your application 
>> code? In 
>> those cases where it crashes, is this_copy_indices.size()==0, so that the 
>> call 
>> to .back() yields nonsense? 
>>
>> As to why this happens, and does so only intermittently, I don't know. A 
>> common cause is that some earlier piece of code (accidentally) overwrites 
>> data. But I can of course only speculate. 
>>
>> 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/66c52a94-54da-4204-9b0d-17d1e4eaa3b8%40googlegroups.com.

Re: [deal.II] Re: MGTransferMatrixFree<> interpolate_to_mg() fails with segfault, when compiled in release mode

['Maxi Miller' via deal.II User Group]

Hei,
the program crashes in release mode way before it crashes in debug mode. 
The code shown is completely from the deal.II-library, with some added 
std::cout<<-lines. The only code written by myself is in main.cc, but the 
crash happens in the code from the library.


Am Donnerstag, 1. August 2019 05:27:21 UTC+2 schrieb Wolfgang Bangerth:
>
>
> Maxi, 
>
> since your program is crashing in debug mode even before the location 
> where it 
> crashes in release mode, let's just ignore the release mode issue -- if a 
> program doesn't run in debug mode, you should not assume that it does 
> anything 
> reasonable at all in release mode. 
>
>
> > When compiling in debug mode, my results vary, and the program behaves 
> > unpredictable, even though I run it in serial mode. In this case it 
> fails 
> > earlier during 
> > 
> > | 
> >  std::cout <<"Looped over all levels, performing plain_copy: 
> > "< > "< > if(perform_plain_copy) 
> > { 
> > // In this case, we can simply copy the local range. 
> > AssertDimension(dst[dst.max_level()].local_size(),src.local_size()); 
> >  dst[dst.max_level()].copy_locally_owned_data_from(src); 
> > return; 
> > } 
> > elseif(perform_renumbered_plain_copy) 
> > { 
> > //std::cout << "dst_level-size: " << 
> dst[dst.max_level()].local_size() 
> > << '\n'; 
> > //std::cout << "src-size: " << src.local_size() << '\n'; 
> > Assert(copy_indices_level_mine.back().empty(),ExcInternalError()); 
> > LinearAlgebra::distributed::Vector_level = 
> >  dst[dst.max_level()]; 
> > for(std::pairi :this_copy_indices.back()) 
> >  
> dst_level.local_element(i.second)=src.local_element(i.first); 
> >  std::cout <<"Renumbered copy done\n"; 
> > return; 
> > } 
> > | 
> > 
> > plain_copy is (at some point) equal to 0, resulting in the program 
> choosing 
> > the second option (perform_renumbered_plain_copy), where it crashes with 
> a 
> > segfault while evaluating the loop 
> > | 
> > 
> > 
> > for(std::pairi :this_copy_indices.back()) 
> >  
> dst_level.local_element(i.second)=src.local_element(i.first); 
> > | 
> > 
> > This behaviour is rather random, and changes at different runs even 
> though the 
> > program was not recompiled. 
> > 
> > Why is that happening? Did I forget to initialize something, or is that 
> an 
> > internal bug? 
>
> Is the code you show somewhere in the library, or in your application 
> code? In 
> those cases where it crashes, is this_copy_indices.size()==0, so that the 
> call 
> to .back() yields nonsense? 
>
> As to why this happens, and does so only intermittently, I don't know. A 
> common cause is that some earlier piece of code (accidentally) overwrites 
> data. But I can of course only speculate. 
>
> 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/0ee8a844-3691-4ad9-b55f-8c50f49863aa%40googlegroups.com.

[deal.II] MGTransferMatrixFree<> interpolate_to_mg() fails with segfault, when compiled in release mode

['Maxi Miller' via deal.II User Group]

I am trying to implement a test program using the matrix-free method, 
combined with a geometric multigrid preconditioner. In my case my 
vmult/apply_add-function needs access to the solution from the last step. 
Thus I added a vector
LinearAlgebra::distributed::Vector old_solution;
which is filled using
template 
void LaplaceOperator::update_old_solution(const 
LinearAlgebra::distributed::Vector ){
this->old_solution = 
LinearAlgebra::distributed::Vector(solution);
}

For each level of grids I need to fill this vector with the current 
solution, interpolated onto the correct mesh. Thus, I added the following 
code:

const unsigned int nlevels = triangulation.n_global_levels();
mg_matrices.resize(0, nlevels - 1);
std::set dirichlet_boundary;
dirichlet_boundary.insert(0);
mg_constrained_dofs.initialize(dof_handler);
mg_constrained_dofs.make_zero_boundary_constraints(dof_handler,
   dirichlet_boundary);

MGTransferMatrixFree mg_transfer(mg_constrained_dofs);

  const unsigned int max_level =
dof_handler.get_triangulation().n_global_levels() - 1;
  const unsigned int min_level = 0;
  MGLevelObject>
level_projection(min_level, max_level);
pcout << "Interpolate to mg\n";
  mg_transfer.interpolate_to_mg(dof_handler, level_projection, solution);
pcout << "Interpolated to mg\n";

for (unsigned int level = 0; level < nlevels; ++level)
{
pcout << "Level: " << level << '\n';
IndexSet relevant_dofs;
DoFTools::extract_locally_relevant_level_dofs(dof_handler,
  level,
  relevant_dofs);
AffineConstraints level_constraints;
level_constraints.reinit(relevant_dofs);
level_constraints.add_lines(
mg_constrained_dofs.get_boundary_indices(level));
level_constraints.close();
typename MatrixFree::AdditionalData additional_data;
additional_data.tasks_parallel_scheme =
MatrixFree::AdditionalData::partition_partition;
additional_data.mapping_update_flags =
(update_gradients | update_JxW_values | 
update_quadrature_points);
additional_data.level_mg_handler = level;
std::shared_ptr> mg_mf_storage_level(
new MatrixFree());
mg_mf_storage_level->reinit(dof_handler,
level_constraints,
QGauss<1>(fe.degree + 1),
additional_data);
mg_matrices[level].initialize(mg_mf_storage_level,
  mg_constrained_dofs,
  level);
mg_matrices[level].set_time_step(time_step);

mg_matrices[level].update_old_solution(level_projection[level]);
}

My program works without crash (even though I get wrong results, but that 
is another problem) when compiling it in debug mode, but crashes 
immediately after pcout << "Interpolate to mg\n" with a segfault; when 
running in release mode. Why is that? Or are there other methods I could 
use for interpolation?
Thanks!

-- 
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/c8e39843-b9d1-4694-ac8c-bef018b42782%40googlegroups.com.
/* -
 *
 * Copyright (C) 2009 - 2019 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, redistribute
 * it, and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * The full text of the license can be found in the file LICENSE.md at
 * the top level directory of deal.II.
 *
 * -
 *
 * Authors: Katharina Kormann, Martin Kronbichler, Uppsala University,
 * 2009-2012, updated to MPI version with parallel vectors in 2016
 */
#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 

#include 
#include 
namespace Step37
{
using namespace dealii;
const unsigned int degree_finite_element = 3;
const unsigned int dimension = 2;

template 
class Coefficient : 

Re: [deal.II] Error on runnig the examples

['Maxi Miller' via deal.II User Group]

Apparently you did not install the library correct, there is no *.so-file 
in that folder. Maybe try to reinstall it?

Am Mittwoch, 31. Juli 2019 09:38:17 UTC+2 schrieb Javad Hashemi:
>
> Dear Dr. Bangreth,
>
> I really appreciate your email. I checked the lib folder and you can see 
> the list of files as follow:
>
> root@Lenovo-PC:~#  ls -l /root/math/dealii-9.1.1/lib/
> total 0
> drwxrwxrwx 1 root root 4096 Jul 28 02:11 cmake
> drwxrwxrwx 1 root root 4096 Jul 28 02:12 pkgconfig
>
> would you please help me to know the reason for the error?
>
> Thanks.
> Sincerely,
>
> Javad 
>
>
>
>
> On Tue, Jul 30, 2019 at 12:19 PM Wolfgang Bangerth  > wrote:
>
>> On 7/30/19 5:41 AM, Javad Hashemi wrote:
>> > All files exist and accessible. 
>>
>> But it doesn't seem to be available when you run cmake. Can you show us 
>> what you get when you do
>>ls -l /root/math/dealii-9.1.1/lib/
>> ? I will note that in general, it's a bad idea to run things as root. 
>> You should be running your programs as a regular user.
>>
>> 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 dea...@googlegroups.com .
>> To view this discussion on the web visit 
>> https://groups.google.com/d/msgid/dealii/4abe20d1-d7a8-4b24-4660-3498d3993a87%40colostate.edu
>> .
>>
>

-- 
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/ca27c645-e19a-469b-86b3-fff73cc980f9%40googlegroups.com.

Re: [deal.II] Access old solution(s) in MatrixFreeOperators::Base::apply_add()

['Maxi Miller' via deal.II User Group]

I found my mistake, I could not reuse the same partition for the main 
matrix and the MG-matrices, thus I had to add 


LinearAlgebra::distributed::Vector local_solution;

mg_matrices[level].initialize_dof_vector(local_solution);

//local_solution = solution;

mg_matrices[level].set_old_solution(local_solution);

in the loop while creating those matrices. My next problem, though, is how 
to fill the local_solution-vector based on the solution-vector. When 
uncommenting the line above, I get exactly the same error as before. Are 
there alternatives?


Am Freitag, 26. Juli 2019 19:36:49 UTC+2 schrieb Daniel Arndt:
>
> Maxi,
>
> [...] I tried adding a class member and initializing it during 
>> setup_system(), but when accessing it via read_dof_values(), I had a layout 
>> mismatch, resulting in
>> The violated condition was: 
>> vec.partitioners_are_compatible(*dof_info.vector_partitioner)
>> Additional information: 
>> The parallel layout of the given vector is not compatible with the 
>> parallel partitioning in MatrixFree. Use MatrixFree::initialize_dof_vector 
>> to get a compatible vector.
>>
>> Did you initialize its layout via MatrixFree::initialize_dof_vector() 
>  in fact?
>
> Best,
> Daniel
>

-- 
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/ae0ba1a3-d2bc-4d95-a5a6-1a89331a0c4f%40googlegroups.com.

[deal.II] Access old solution(s) in MatrixFreeOperators::Base::apply_add()

['Maxi Miller' via deal.II User Group]

I followed step-37 for creating a MatrixFreeOperator, but now I also would 
like to access a previous solution (which is needed during apply_add()), 
and stored using a LinearAlgebra::distributed::Vector. In step-48 there is 
an example which uses previous solutions, but here the solve-function is 
replaced, and therefore making it difficult to implement it in my problem.
Thus, what would be the most efficient way of giving the operator-function 
apply_add access to those solution vectors? I tried adding a class member 
and initializing it during setup_system(), but when accessing it via 
read_dof_values(), I had a layout mismatch, resulting in
The violated condition was: 
vec.partitioners_are_compatible(*dof_info.vector_partitioner)
Additional information: 
The parallel layout of the given vector is not compatible with the 
parallel partitioning in MatrixFree. Use MatrixFree::initialize_dof_vector 
to get a compatible vector.

Are there better ways?
Thanks!

-- 
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/855a446f-4bde-4a56-8ae2-1871488f0b65%40googlegroups.com.

Re: [deal.II] VMult-function in the MatrixFree-Framework

['Maxi Miller' via deal.II User Group]

Thus I intended to use a GMG-preconditioner, with a Chebyshev-smoother. 
Usually the GMG-preconditioner does not require the main matrix either 
(yes, small matrices are still required, but they do not require the same 
amount of storage space), but here I run into the problem which is 
mentioned in a separate question, i.e. I can not get the gradients/values 
of my function on non-active cells while having to loop over both 
non-active and active cells for building the MG-matrices, if that makes 
sense?

Am Mittwoch, 24. Juli 2019 21:18:13 UTC+2 schrieb Wolfgang Bangerth:
>
> On 7/24/19 12:23 PM, 'Maxi Miller' via deal.II User Group wrote: 
> > My LinearOperator only provides a vmult-interface, nothing else, but for 
> > initializing the preconditioner I still need a sparse matrix, thus I can 
> not 
> > use the LinearOperator, as far as I understand. Or is there a way to use 
> it? 
>
> You just found the problem everyone faces. The LinearOperator interface is 
> useful if you only know how to express a linear operator through its 
> action, 
> without knowing the entries of the matrix that represent this action. This 
> is 
> all you need for the iterative solvers we have. 
>
> But to work well, these iterators all need preconditioners, and most of 
> the 
> preconditioners we have require knowledge of the entries of the matrix -- 
> i.e., there is no longer any advantage to using LinearOperator over 
> building a 
> matrix. The ways you can work around this are Chebyshev iterations, for 
> example; step-20 also in great detail explores this issue. But the point 
> is 
> that it's a difficult issue. 
>
> 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/6f2c93b2-f316-4a99-b86b-4775d450b2dd%40googlegroups.com.

Re: [deal.II] VMult-function in the MatrixFree-Framework

['Maxi Miller' via deal.II User Group]

My LinearOperator only provides a vmult-interface, nothing else, but for 
initializing the preconditioner I still need a sparse matrix, thus I can 
not use the LinearOperator, as far as I understand. Or is there a way to 
use it?

Am Mittwoch, 24. Juli 2019 17:28:14 UTC+2 schrieb Daniel Arndt:
>
> Is there then a way to use LinearOperator as Input for a preconditioner? 
>> Else I still have to form the full system matrix (which I would like to 
>> avoid).
>>
>
> Most precoditioners require access to individual matrx entries. If the 
> class you derive from linear operator provides a suitable interface, you 
> can use it in a preconditioner.
> Building suitable preconditioners compatible with matrix-free methods is 
> non-trivial. The default choice would be PreconditionChebyshev.
>
> Best,
> Daniel
>

-- 
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/34ca9d5a-5a0c-408a-8479-8d17fcb653f2%40googlegroups.com.

Re: [deal.II] Matrix-free-method for an operator depending on input from solver

['Maxi Miller' via deal.II User Group]

Hei,
based on what I know from literature a GMG-preconditioner with Chebyshev 
smoothing should work rather well for my case, thus I intended to follow 
example 37 here. If I can set v_i to 1 for calculating the diagonal, that 
would simplify the problem quite a lot.
The main reason for testing the matrix-free approach in my case is enormous 
RAM-consumption when creating the matrix (I have to use a rather dense 
grid, else there will be errors), and a long build time for the matrix 
itself.
I hope that makes sense?

Am Mittwoch, 24. Juli 2019 17:51:48 UTC+2 schrieb Martin Kronbichler:
>
> Dear Maxi,
>
> There is not really a simple answer to your request: It depends on how 
> good a preconditioner you need. Matrix-free methods work best if you can 
> use simple preconditioners in the sense that they only need some matrix 
> entries (if any) and the matrix-vector product (or operator evaluation in 
> the nonlinear case). Note that I also consider multigrid with Chebyshev 
> smoothing as in step-37 and step-59 a simple preconditioner. It is up to 
> you to find out whether your problem admits such a simple form and 
> iterative solvers, such that GMRES, some nonlinear Newton-Krylov or 
> nonlinear Krylov variant, can cope with that and not use too many 
> iterations. We as a mailing list cannot help with the preconditioner 
> because this is a research question that many people would like an answer 
> of. My experience is that matrix-free setups are mostly done for 
> performance reasons (because the mat-vec is cheaper) for well-understood 
> systems, like the Laplacian with multigrid or some linear/nonlinear 
> time-dependent equations with dominating mass matrix terms where simple 
> iterative methods suffice. Sometimes matrix-free is done because one can 
> fit larger problems in memory. But all cases where I know it is used are 
> done when one at least has an expectation of what the matrix and its 
> eigenvalue spectrum look like. If you do not have that, because you have 
> some difficult linear/nonlinear systems, matrix-free setups might help 
> provide you with a black-box view - but they make the preconditioner 
> selection harder, not easier, because you restrict yourself in the choice 
> of preconditioner to the more basic variants.
>
> Now to the technical part: If you need the diagonal of a matrix and rely 
> on an accurate representation, you will need to build the linearization in 
> some way or the other. This is also what step-37 or step-59 do - they 
> compute the local stiffness matrix, throw away everything but the diagonal, 
> and use that inside some solver. Tor the diagonal you would run through all 
> unit vectors with v, i.e., you add (1,0,0,...), (0,1,0,...), and so on. 
> Then you only keep the one entry where v was one. Your approach to do this 
> globally is certainly not going to scale because it is an O(n^2) operation 
> to check every unit vector. In principle, you could do it locally 
> element-by-element similar to step-37 also for nonlinear operators. The 
> question is whether you have the full action inside a single "local_apply" 
> function such that you could start adding unit vectors there. Given a 
> nonlinear field u that you linearize about, you could perturb it with e*v_i 
> and compute the local diagonal according to your formula. From what I 
> understand, you have the operator in an even more abstract way, so 
> collecting everything in one local call with matrix-free facilities would 
> also be non-trivial.
>
> Best,
> Martin
> On 24.07.19 15:01, 'Maxi Miller' via deal.II User Group wrote:
>
> I am currently trying to implement the JFNK-method in the matrix-free 
> framework (by following step-37 and step-59) for solving nonlinear 
> equations of the form F(u)=0. The method itself replaces the multiplication 
> of the explicit jacobian J with the krylov vector v during solving with the 
> approximation (F(u+e*v)-F(u))/e, which removes the matrix. 
> The initial test using a LinearOperator worked, where I modified the 
> vmult-function accordingly (c.f. my earlier questions here), but still I 
> had to form the jacobian at least once (which I wanted to avoid), thus the 
> test of the matrix-free framework. 
> In theory I just should have to provide the vmult-function, but based on 
> step-37 and step-59 I also have to provide the full operator at least once 
> during initialization, as f.ex. to form the inverse diagonal in step-37. 
> This would be the same as calculating (F(u_i+e*v_i)-F(u_i))/e, but I do not 
> know the value of v_i.
> Thus, what would be the best approach of implementing this operator using 
> the MatrixFree-framework in deal.II, without having to form a single 
> matrix, and without knowing the exact value of F? Ideal

[deal.II] Matrix-free-method for an operator depending on input from solver

['Maxi Miller' via deal.II User Group]

I am currently trying to implement the JFNK-method in the matrix-free 
framework (by following step-37 and step-59) for solving nonlinear 
equations of the form F(u)=0. The method itself replaces the multiplication 
of the explicit jacobian J with the krylov vector v during solving with the 
approximation (F(u+e*v)-F(u))/e, which removes the matrix. 
The initial test using a LinearOperator worked, where I modified the 
vmult-function accordingly (c.f. my earlier questions here), but still I 
had to form the jacobian at least once (which I wanted to avoid), thus the 
test of the matrix-free framework. 
In theory I just should have to provide the vmult-function, but based on 
step-37 and step-59 I also have to provide the full operator at least once 
during initialization, as f.ex. to form the inverse diagonal in step-37. 
This would be the same as calculating (F(u_i+e*v_i)-F(u_i))/e, but I do not 
know the value of v_i.
Thus, what would be the best approach of implementing this operator using 
the MatrixFree-framework in deal.II, without having to form a single 
matrix, and without knowing the exact value of F? Ideally it can be used 
for a preconditioner too, such as the GMG in step-37 and step-59.
Thanks!

-- 
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/257bb96a-ba19-48cf-a7c4-49e6f178e417%40googlegroups.com.

Re: [deal.II] VMult-function in the MatrixFree-Framework

['Maxi Miller' via deal.II User Group]

Is there then a way to use LinearOperator as Input for a preconditioner? 
Else I still have to form the full system matrix (which I would like to 
avoid).

Am Mittwoch, 24. Juli 2019 13:43:09 UTC+2 schrieb Daniel Arndt:
>
> On the other hand, would it be sufficient to unroll the 
 residual-calculation, i.e. if the residual is F(u) = nabla^2 u + f, to 
 write (in the cell-loop): (nabla^2(u + eps*src) + f - nabla^2u - 
 f)/epsilon? Or would that lead to wrong results?

>>>  
>>> It is sufficient to only change compute_diagonal() to use the MatrixFree 
>>> framework. Passing multiple vectors at once might be a little bit more 
>>> difficult to implement.
>>>
>> Could you elaborate on that further (changing the diagonal value)?  
>>
>
> Oh, that was a typo. I meant that you only have to implement 
> calculate_residual and can leave the rest of you LinearOperator as is.
>
> Best,
> Daniel
>

-- 
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/faa5a96a-683c-4c90-a3a9-3f1d0bd38be1%40googlegroups.com.

Re: [deal.II] VMult-function in the MatrixFree-Framework

['Maxi Miller' via deal.II User Group]

Am Mittwoch, 24. Juli 2019 00:39:14 UTC+2 schrieb Daniel Arndt:
>
> Maxi,
>
> No, that code is the vmult-call I am using in my LinearOperator (and is 
>> used directly in my solve()-call). The calculate_residual()-function is 
>> similar to the function in step-15, with the difference, that I do not 
>> provide alpha, but the vector, and return the calculated residual value.
>>
>  
> So this is where your cell loop currently lives and that can be 
> represented by a matrix-vector product. :-) 
>  
>
>> Thus, as far as I could understand, I have to do the following in the 
>> matrix-free cell loop:
>>
>>- Calculate residual of the current solution
>>- Calculate residual of the current solution + epsilon * src
>>- Calculate dst, and return it
>>
>> On the other hand, would it be sufficient to unroll the 
>> residual-calculation, i.e. if the residual is F(u) = nabla^2 u + f, to 
>> write (in the cell-loop): (nabla^2(u + eps*src) + f - nabla^2u - 
>> f)/epsilon? Or would that lead to wrong results?
>>
>  
> It is sufficient to only change compute_diagonal() to use the MatrixFree 
> framework. Passing multiple vectors at once might be a little bit more 
> difficult to implement.
>
Could you elaborate on that further (changing the diagonal value)?  

>
> Currently, you residual seems to depend linearly on u. In that case, the 
> difference does not depend on u at all.
>
I know, the equation is linear anyway, but was intended to be used as 
simple example. For further tests the non-linear equations will be used.

>
> Best,
> Daniel
>

-- 
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/12247d99-b2a8-4125-b98b-6e965e9f6f4d%40googlegroups.com.

Re: [deal.II] Re: Matrix-Vector-multiplication on deal.II-grid

['Maxi Miller' via deal.II User Group]

I took a look at your (incomplete) step-58, which (as far as I could 
understand) solves the equation, but without propagating it in space, just 
in time. I would like to propagate the result in time and accordingly in 
space, i.e. if the pulse propagates through a fiber, I would like to see 
the moving pulse. Is that possible?
The matrix I was mentioning earlier would be a dense matrix, formed 
according to the DHT-method (i.e. no relation to FEM)

Am Dienstag, 23. Juli 2019 22:35:57 UTC+2 schrieb Wolfgang Bangerth:
>
> On 7/22/19 1:32 AM, 'Maxi Miller' via deal.II User Group wrote: 
> > More or less, yes, together with how I could get the "column" vector I 
> would 
> > like to apply the matrix to. It is not a problem to do that in a 
> FTDT-setting, 
> > but I do not know how to do that in a FEM-setting. 
>
> When you say "column vector", you mean the vector of unknowns that 
> corresponds 
> to one column of nodes in the mesh? 
>
Yes, basically 

>
> One can think of deal.II Triangulations as unstructured meshes. So 
> extracting 
> a column of nodes is not a natural thing to do. But you can of course 
> always 
> loop over all nodes (e.g., using DoFTools::map_dofs_to_support_points()) 
> and 
> group them by the x-coordinate of their location. 
>
> How you would build this matrix? That's something we can't know without 
> actually seeing a formula for the operator you're thinking of. 
>
> 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/4497b9c1-b6f9-4d51-b917-762a13004423%40googlegroups.com.

Re: [deal.II] VMult-function in the MatrixFree-Framework

['Maxi Miller' via deal.II User Group]

No, that code is the vmult-call I am using in my LinearOperator (and is 
used directly in my solve()-call). The calculate_residual()-function is 
similar to the function in step-15, with the difference, that I do not 
provide alpha, but the vector, and return the calculated residual value. 
Thus, as far as I could understand, I have to do the following in the 
matrix-free cell loop:

   - Calculate residual of the current solution
   - Calculate residual of the current solution + epsilon * src
   - Calculate dst, and return it

On the other hand, would it be sufficient to unroll the 
residual-calculation, i.e. if the residual is F(u) = nabla^2 u + f, to 
write (in the cell-loop): (nabla^2(u + eps*src) + f - nabla^2u - 
f)/epsilon? Or would that lead to wrong results?

Am Dienstag, 23. Juli 2019 22:31:06 UTC+2 schrieb Daniel Arndt:
>
> My aim is to implement the jacobian approximation, i.e. (if F(u) = 0 is 
>> the function to solve) I wanted to implement
>> dst = (F(current_solution + epsilon * src) - F(current_solution))/epsilon.
>> In my LinearOperator, my code is 
>> LinearOperator jacobian_operator;
>> jacobian_operator.vmult = [&](LinearAlgebraTrilinos::MPI::Vector , 
>> const LinearAlgebraTrilinos::MPI::Vector ){
>> local_src = src;
>> extended_src = src;
>> local_solution = present_solution;
>>
>> double epsilon = 1e-6;
>> if(local_src.l2_norm() != 0){
>> epsilon = sqrt((1 + local_solution.l2_norm()) * 1e-16) / 
>> local_src.l2_norm();
>> };
>> forward_eps_solution = present_solution;
>> backward_eps_solution = present_solution;
>>
>> extended_src *= epsilon;
>> forward_eps_solution += extended_src;
>> backward_eps_solution -= extended_src;
>>
>> compute_residual(backward_eps_solution, cur_residual);
>> compute_residual(forward_eps_solution, eps_residual);
>> eps_residual -= cur_residual;
>> eps_residual /= (2 * epsilon);
>> dst.reinit(locally_owned_dofs,
>>mpi_communicator);
>> dst = eps_residual;
>> };
>>
>> Is the same/similar possible for the matrix-free-approach?
>>
> Sure, you can always write a wrapper around the actual vmult call (which I 
> assume is happening in compute_residual()).
>
> Best,
> Daniel
>
>

-- 
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/5f0807fd-6829-4d60-b49f-41b7e1bb2e27%40googlegroups.com.

Re: [deal.II] VMult-function in the MatrixFree-Framework

['Maxi Miller' via deal.II User Group]

My aim is to implement the jacobian approximation, i.e. (if F(u) = 0 is the 
function to solve) I wanted to implement
dst = (F(current_solution + epsilon * src) - F(current_solution))/epsilon.
In my LinearOperator, my code is 
LinearOperator jacobian_operator;
jacobian_operator.vmult = [&](LinearAlgebraTrilinos::MPI::Vector , const 
LinearAlgebraTrilinos::MPI::Vector ){
local_src = src;
extended_src = src;
local_solution = present_solution;

double epsilon = 1e-6;
if(local_src.l2_norm() != 0){
epsilon = sqrt((1 + local_solution.l2_norm()) * 1e-16) / local_src.
l2_norm();
};
forward_eps_solution = present_solution;
backward_eps_solution = present_solution;

extended_src *= epsilon;
forward_eps_solution += extended_src;
backward_eps_solution -= extended_src;

compute_residual(backward_eps_solution, cur_residual);
compute_residual(forward_eps_solution, eps_residual);
eps_residual -= cur_residual;
eps_residual /= (2 * epsilon);
dst.reinit(locally_owned_dofs,
   mpi_communicator);
dst = eps_residual;
};

Is the same/similar possible for the matrix-free-approach?
Thanks!

Am Dienstag, 23. Juli 2019 19:00:14 UTC+2 schrieb Daniel Arndt:
>
>
>
>>- When using the LinearOperator, the function vmult() is used when 
>>solving the system with a GMRES-solver (or similar). Does the same hold 
>> up 
>>for the MatrixFree-framework?
>>
>> Yes. 
>
>>
>>- If the above is correct: In each GMRES-iteration I have to 
>>calculate the residual of my equation, once for the current solution and 
>>once for the current solution + a small step. At the moment I calculate 
>> the 
>>residual by summing the terms into my residual vector, and then 
>>distributing the values while setting the boundary condition to 0. If I 
>>understand it correctly, that happens as soon as I call 
>>phi.distribute_local_to_global(dst), with phi a FEEvaluation-element.
>>
>> Yes. 
>
>>
>>- This means that I first have to assemble the residual vector, then 
>>distribute it, and then can use the vmult-function. Is that correct, or 
>> do 
>>I have to use another approach?
>>
>> Isn't assembling and distributing already what `vmult` is supposed to do? 
> Can you elaborate on what you are doing in each GMRES iteration?
>
> Best,
> Daniel
>
>

-- 
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/0240a020-581e-4202-bdc4-57ca8535a182%40googlegroups.com.

[deal.II] VMult-function in the MatrixFree-Framework

['Maxi Miller' via deal.II User Group]

In my project I already managed to replace the system_matrix in the 
solve()-call by a custom LinearOperator, but still I had to retain the 
matrix itself (for the preconditioner, for example). Thus I am now trying 
to use the MatrixFree-framework for that.
Thus, some questions arose:

   - When using the LinearOperator, the function vmult() is used when 
   solving the system with a GMRES-solver (or similar). Does the same hold up 
   for the MatrixFree-framework?
   - If the above is correct: In each GMRES-iteration I have to calculate 
   the residual of my equation, once for the current solution and once for the 
   current solution + a small step. At the moment I calculate the residual by 
   summing the terms into my residual vector, and then distributing the values 
   while setting the boundary condition to 0. If I understand it correctly, 
   that happens as soon as I call phi.distribute_local_to_global(dst), with 
   phi a FEEvaluation-element. This means that I first have to assemble the 
   residual vector, then distribute it, and then can use the vmult-function. 
   Is that correct, or do I have to use another approach?

Thanks!

-- 
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/a2dd198f-513d-4ee0-9ebc-2da339b4f133%40googlegroups.com.

Re: [deal.II] Code tries to access DoF indices on artificial cells, even though those cells should not be artificial

['Maxi Miller' via deal.II User Group]

Furthermore, based on my understanding, I have to loop over all cells in 
order to assemble the multigrid-matrices, but can only aquire the gradients 
on the current cell if it is active. Wouldn't that crash (looping and 
gathering data from all cells, but only getting data from the active cells)?

Am Samstag, 20. Juli 2019 11:03:05 UTC+2 schrieb Maxi Miller:
>
> Based on what I get (within the if-loop):
> I print the cell information using 
> std::cout << "Aquiring function gradients on cell with data points " << 
> cell->level_subdomain_id() << '\t' << triangulation.
> locally_owned_subdomain() << '\t' << cell->active() << '\t' << 
> cell_counter << '\n';
>
> with cell_counter a static variable, increasing for every function call.
> This gives me
> Aquiring function gradients on cell with data points 1  1   0   1
> Aquiring function gradients on cell with data points 0  0   0   1
> Aquired function gradients on cell 1
> Aquired function gradients on cell 1
> Aquiring function gradients on cell with data points 1  1   0   2
> Aquiring function gradients on cell with data points 0  0   0   2
> Aquired function gradients on cell 2
> Aquired function gradients on cell 2
> Aquiring function gradients on cell with data points 0  0   0   3
> Aquiring function gradients on cell with data points 1  1   0   3
> Aquired function gradients on cell 3
>
>
> before crashing, indicating a problem on cell 1 1 0 3. As long as 
> cell->active() is false, I can not call the function is_artificial() for 
> checking if I am on an artificial cell. Is my if-loop correct, or did I 
> miss something here?
>
> Am Freitag, 19. Juli 2019 18:01:03 UTC+2 schrieb Wolfgang Bangerth:
>>
>> On 7/19/19 9:59 AM, 'Maxi Miller' via deal.II User Group wrote: 
>> > What do the parallel multigrid functions do? Could you just have 
>> your if() 
>> > inside that function -- i.e., it is called on all cells, but does 
>> nothing if 
>> > the cell doesn't allow you to do anything? 
>> > 
>> > 
>> > That's exactly what my if()-conditions should do (but is not doing)... 
>>
>> Then put a breakpoint on it and see why it isn't doing what you think it 
>> is 
>> doing :-) 
>>
>> 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/27a58ff4-e4fb-4510-8220-68d4979f346a%40googlegroups.com.

Re: [deal.II] Re: Matrix-Vector-multiplication on deal.II-grid

['Maxi Miller' via deal.II User Group]

More or less, yes, together with how I could get the "column" vector I 
would like to apply the matrix to. It is not a problem to do that in a 
FTDT-setting, but I do not know how to do that in a FEM-setting.

Am Montag, 22. Juli 2019 00:05:45 UTC+2 schrieb Wolfgang Bangerth:
>
> On 7/12/19 11:28 AM, 'Maxi Miller' via deal.II User Group wrote: 
> > It is difficult to write it as a single integral. The operation is 
> similar to 
> > the split-step fourier method, i.e. transforming the column vector f(r) 
> once using 
> > g(rho)=2\pi\int_0^\infty rf(r)J_0(2\pi\rho r)dr, 
> > multiplying it with a vector, and transforming it back using 
> > f(r) = 2\pi\int_0^\infty\rho g(\rho)J_0(2\pi\rho r)d\rho 
> > The operation is for radially symmetric systems, i.e. with z along the 
> x-axis, 
> > and r along the y-axis. When starting on the left border with f_0, i.e. 
> at 
> > position z = 0, doing the operation mentioned above gives the values for 
> the 
> > nodes at z = 1, when enumerating the nodes from 0 to n along the z axis, 
> and 
> > having equidistant nodes along the z-axis. Those integrals can be 
> replaced by 
> > a matrix-vector-multiplication, thus making it easier to implement 
> numerically. 
>
> This makes sense -- every linear operator can be represented as a matrix. 
> So 
> yes, you can express the operation as a matrix, and of course you can 
> express 
> this matrix using the deal.II classes. 
>
> But I'm not clear what your question is then. Are you asking how to build 
> this 
> matrix? 
>
> 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/c11c3cf6-4f75-4fb4-abfb-85f5f61c28ba%40googlegroups.com.

Re: [deal.II] Code tries to access DoF indices on artificial cells, even though those cells should not be artificial

['Maxi Miller' via deal.II User Group]

Based on what I get (within the if-loop):
I print the cell information using 
std::cout << "Aquiring function gradients on cell with data points " << cell
->level_subdomain_id() << '\t' << triangulation.locally_owned_subdomain() << 
'\t' << cell->active() << '\t' << cell_counter << '\n';

with cell_counter a static variable, increasing for every function call.
This gives me
Aquiring function gradients on cell with data points 1  1   0   1
Aquiring function gradients on cell with data points 0  0   0   1
Aquired function gradients on cell 1
Aquired function gradients on cell 1
Aquiring function gradients on cell with data points 1  1   0   2
Aquiring function gradients on cell with data points 0  0   0   2
Aquired function gradients on cell 2
Aquired function gradients on cell 2
Aquiring function gradients on cell with data points 0  0   0   3
Aquiring function gradients on cell with data points 1  1   0   3
Aquired function gradients on cell 3


before crashing, indicating a problem on cell 1 1 0 3. As long as 
cell->active() is false, I can not call the function is_artificial() for 
checking if I am on an artificial cell. Is my if-loop correct, or did I 
miss something here?

Am Freitag, 19. Juli 2019 18:01:03 UTC+2 schrieb Wolfgang Bangerth:
>
> On 7/19/19 9:59 AM, 'Maxi Miller' via deal.II User Group wrote: 
> > What do the parallel multigrid functions do? Could you just have 
> your if() 
> > inside that function -- i.e., it is called on all cells, but does 
> nothing if 
> > the cell doesn't allow you to do anything? 
> > 
> > 
> > That's exactly what my if()-conditions should do (but is not doing)... 
>
> Then put a breakpoint on it and see why it isn't doing what you think it 
> is 
> doing :-) 
>
> 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/466b7db8-a8f1-4215-a5c1-41ba16132e76%40googlegroups.com.

Re: [deal.II] Create a LinearOperator with a function instead of a matrix

['Maxi Miller' via deal.II User Group]

Would that work together with the GMG-preconditioner? Will test it there as 
soon as my problem in my second question could be solved.
There I have the following code:

SolverControlcoarse_solver_control(1000, 1e-10, false, false);
SolverGMRES   coarse_solver(coarse_solver_control);
PreconditionIdentity id;
MGCoarseGridIterativeSolver,
matrix_t,
PreconditionIdentity>
coarse_grid_solver(coarse_solver, coarse_matrix, id);

using Smoother = LA::MPI::PreconditionJacobi;
MGSmootherPrecondition mg_smoother;
mg_smoother.initialize(mg_matrices, Smoother::AdditionalData(0.5));
mg_smoother.set_steps(2);

mg::Matrix mg_matrix(mg_matrices);
mg::Matrix mg_interface_up(mg_interface_matrices);
mg::Matrix mg_interface_down(mg_interface_matrices);

Multigrid mg(mg_matrix, coarse_grid_solver, mg_transfer, 
mg_smoother, mg_smoother);
mg.set_edge_matrices(mg_interface_down, mg_interface_up);

PreconditionMG> 
preconditioner(dof_handler,

   mg,

   mg_transfer);

and here I intent to replace all occurences of matrix_t with my custom 
LinearOperator (if that makes sense).

Am Freitag, 19. Juli 2019 17:51:10 UTC+2 schrieb Wolfgang Bangerth:
>
> On 7/19/19 9:47 AM, 'Maxi Miller' via deal.II User Group wrote: 
> > What I did was the replacement of system_matrix with a function, 
> providing a 
> > vmult()-function. Now I would like to add a preconditioner based on that 
> > function to reduce the necessary GMRES-iterations. But until now most of 
> the 
> > preconditioners require a matrix for initialization. Thus I wanted to 
> replace 
> > that matrix with my LinearOperator, if that is possible? 
>
> Yes, as I already said, this is possible. Have you tried? 
>
> I think you should also be able to look at step-20 (and maybe 22) for 
> examples. 
>
> 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/0c8f71f7-8d29-4a4f-abb1-65afbcc72669%40googlegroups.com.

Re: [deal.II] Code tries to access DoF indices on artificial cells, even though those cells should not be artificial

['Maxi Miller' via deal.II User Group]

>
> What do the parallel multigrid functions do? Could you just have your if() 
> inside that function -- i.e., it is called on all cells, but does nothing 
> if 
> the cell doesn't allow you to do anything? 
>

That's exactly what my if()-conditions should do (but is not doing)...
Am Freitag, 19. Juli 2019 17:53:05 UTC+2 schrieb Wolfgang Bangerth:
>
> On 7/19/19 9:49 AM, 'Maxi Miller' via deal.II User Group wrote: 
> > I can not use that function, else I will get 
> > An error occurred in line <3543> of file 
> >  in function 
> >  bool dealii::CellAccessor::is_locally_owned() const 
> [with 
> > int dim = 2; int spacedim = 2] 
> > The violated condition was: 
> >  this->active() 
> > Additional information: 
> >  is_locally_owned() can only be called on active cells! 
>
> Ah, I see. I missed that part. 
>
>
> > In my multigrid-assembly-function I have to loop over all cells, not 
> only the 
> > currently active cells (as it is stated in example 50: /we don't just 
> loop 
> > over all active cells, but in fact all cells, active or not/). Thus I 
> end up 
> > in that function. 
>
> Hm, that's a question for the parallel multigrid folks then. 
>
> What do the parallel multigrid functions do? Could you just have your if() 
> inside that function -- i.e., it is called on all cells, but does nothing 
> if 
> the cell doesn't allow you to do anything? 
>
> 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/3a65f070-a3b5-479b-acf7-f41272771b0f%40googlegroups.com.

Re: [deal.II] Code tries to access DoF indices on artificial cells, even though those cells should not be artificial

['Maxi Miller' via deal.II User Group]

I can not use that function, else I will get 
An error occurred in line <3543> of file 
 in function
bool dealii::CellAccessor::is_locally_owned() const 
[with int dim = 2; int spacedim = 2]
The violated condition was: 
this->active()
Additional information: 
is_locally_owned() can only be called on active cells!

In my multigrid-assembly-function I have to loop over all cells, not only 
the currently active cells (as it is stated in example 50: *we don't just 
loop over all active cells, but in fact all cells, active or not*). Thus I 
end up in that function.

Am Freitag, 19. Juli 2019 17:37:01 UTC+2 schrieb Wolfgang Bangerth:
>
> On 7/19/19 9:31 AM, 'Maxi Miller' via deal.II User Group wrote: 
> > Based on the stacktrace the problematic access happens within 
> mg_cell_worker(): 
>  > ... 
>  > which is guarded by the if-condition. 
>
> Then we have a miracle :-) You'll have to check in a debugger how you got 
> into 
> this place. If you go up to frame #8 (assemble_multigrid) you should be 
> able 
> why the function mg_cell_worker() was called and what the condition is 
> that 
> you checked/wanted to check. 
>
> BTW, I typically use 
>cell->is_locally_owned() 
> in place of the more complicated condition you use. 
>
> 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/d78faba2-5bf4-4c49-a40b-a46eb7d1eb76%40googlegroups.com.

Re: [deal.II] Create a LinearOperator with a function instead of a matrix

['Maxi Miller' via deal.II User Group]

What I did was the replacement of system_matrix with a function, providing 
a vmult()-function. Now I would like to add a preconditioner based on that 
function to reduce the necessary GMRES-iterations. But until now most of 
the preconditioners require a matrix for initialization. Thus I wanted to 
replace that matrix with my LinearOperator, if that is possible?

-- 
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/953da0b8-eb40-4ac1-8316-042c382f2978%40googlegroups.com.

Re: [deal.II] Create a LinearOperator with a function instead of a matrix

['Maxi Miller' via deal.II User Group]

I.e. (based on step-37) I can write something like
LinearOperator jacobian_operator;
//vmult-operation is added later

 using SystemMatrixType = jacobian_operator;
 SystemMatrixType system_matrix;

and follow the example accordingly?

Am Freitag, 19. Juli 2019 17:26:31 UTC+2 schrieb Wolfgang Bangerth:
>
> On 7/19/19 9:18 AM, 'Maxi Miller' via deal.II User Group wrote: 
> > Is it then possible to use this LinearOperator in a preconditioner? The 
> code 
> > works without one (i.e. PreconditionIdentity()), but the Iterations are 
> > increasing quite fast, thus negating all speedup gained from the 
> LinearOperator. 
> > I tried to find something in the examples, but only found data about 
> using 
> > matrices as input for preconditioners, not pure LinearOperators 
> containing 
> > only a vmult-function. 
>
> Yes, this is supposed to work! 
> 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/a77e57c0-c126-4b81-a41b-f95b45342897%40googlegroups.com.

Re: [deal.II] Code tries to access DoF indices on artificial cells, even though those cells should not be artificial

['Maxi Miller' via deal.II User Group]

CopyData&)#1}, 
std::function, 
dealii::TriaIterator, 
true> >, Step15::ScratchData<2>, 
Step15::CopyData>(dealii::TriaIterator, true> > const&, 
dealii::identity, true> > >::type const&, 
dealii::MeshWorker::mesh_loop, true> >, Step15::ScratchData<2>, Step15::CopyData, 
dealii::TriaIterator, 
true> > 
>(dealii::TriaIterator, 
true> > const&, 
dealii::identity, true> > >::type const&, dealii::identity, 
true> > const&, Step15::ScratchData<2>&, Step15::CopyData&)> >::type 
const&, dealii::identity 
>::type const&, Step15::ScratchData<2> const&, Step15::CopyData const&, 
dealii::MeshWorker::AssembleFlags, dealii::identity, 
true> > const&, unsigned int, Step15::ScratchData<2>&, Step15::CopyData&)> 
>::type const&, dealii::identity, 
true> > const&, unsigned int, unsigned int, 
dealii::TriaIterator, 
true> > const&, unsigned int, unsigned int, Step15::ScratchData<2>&, 
Step15::CopyData&)> >::type const&, unsigned int, unsigned 
int)::{lambda(dealii::TriaIterator, true> > const&, Step15::ScratchData<2>&, Step15::CopyData&)#1}, 
std::function, Step15::ScratchData<2> 
const&, Step15::CopyData const&, unsigned int, unsigned int)
#13  ./MG_test: void 
dealii::MeshWorker::mesh_loop, true> >, Step15::ScratchData<2>, Step15::CopyData, 
dealii::TriaIterator, 
true> > 
>(dealii::TriaIterator, 
true> > const&, 
dealii::identity, true> > >::type const&, dealii::identity, 
true> > const&, Step15::ScratchData<2>&, Step15::CopyData&)> >::type 
const&, dealii::identity 
>::type const&, Step15::ScratchData<2> const&, Step15::CopyData const&, 
dealii::MeshWorker::AssembleFlags, dealii::identity, 
true> > const&, unsigned int, Step15::ScratchData<2>&, Step15::CopyData&)> 
>::type const&, dealii::identity, 
true> > const&, unsigned int, unsigned int, 
dealii::TriaIterator, 
true> > const&, unsigned int, unsigned int, Step15::ScratchData<2>&, 
Step15::CopyData&)> >::type const&, unsigned int, unsigned int)
#14  ./MG_test: Step15::MinimalSurfaceProblem<2>::assemble_multigrid()
#15  ./MG_test: Step15::MinimalSurfaceProblem<2>::run()
#16  ./MG_test: main

which is guarded by the if-condition.


Am Freitag, 19. Juli 2019 17:26:00 UTC+2 schrieb Wolfgang Bangerth:
>
> On 7/19/19 9:15 AM, 'Maxi Miller' via deal.II User Group wrote: 
> > I am trying to port example 15 to a MPI-supported version, using the 
> geometric 
> > multigrid method as preconditioner (to get more familiar with the 
> > implementation). I followed example 50 while doing that, but currently I 
> get 
> > the error 
> > An error occurred in line <3866> of file 
> > 
> <~/Downloads/git-files/dealii/include/deal.II/dofs/dof_accessor.templates.h> 
>
> > in function 
> >  void dealii::DoFCellAccessor lda>::get_dof_values(const 
> > InputVector&, ForwardIterator, ForwardIterator) const [with InputVector 
> = 
> > dealii::TrilinosWrappers::MPI::Vector; ForwardIterator = double*; 
> > DoFHandlerType = dealii::DoFHandler<2, 2>; bool level_dof_access = true] 
> > The violated condition was: 
> >  this->is_artificial() == false 
> > Additional information: 
> >  Can't ask for DoF indices on artificial cells. 
> > as soon as I try to run my program using several MPI threads. I tried to 
> guard 
> > that part using if (cell->level_subdomain_id() == 
> > triangulation.locally_owned_subdomain()) as suggested in step-50, but 
> that did 
> > not help. Why am I still trying to request data from artificial cells, 
> even 
> > though I guarded that code part (at least I assume it guards it 
> accordingly). 
>
> Are you sure you guarded the right place? Have you checked with a debugger 
> where the problem actually happens? 
>
> Best 
>   WB 
>
> -- 
>  
> 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/7ec5fd7e-5788-4627-b5b2-4a7bca380981%40googlegroups.com.

Re: [deal.II] Create a LinearOperator with a function instead of a matrix

['Maxi Miller' via deal.II User Group]

Is it then possible to use this LinearOperator in a preconditioner? The 
code works without one (i.e. PreconditionIdentity()), but the Iterations 
are increasing quite fast, thus negating all speedup gained from the 
LinearOperator.
I tried to find something in the examples, but only found data about using 
matrices as input for preconditioners, not pure LinearOperators containing 
only a vmult-function.

Am Dienstag, 16. Juli 2019 18:21:28 UTC+2 schrieb Matthias Maier:
>
> A short comment: 
>
> On Tue, Jul 16, 2019, at 11:18 CDT, Matthias Maier  > wrote: 
>
> > struct LeftVector { 
> > }; 
>
> > struct RightVector { 
> > }; 
>
> These two classes are of course just decoration (showing the minimal 
> interface a vector has to possess). There is usually no need to define 
> custom Vector clases. 
>

-- 
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/c7bc4fb5-1cf7-46ce-af94-2bd5a796715d%40googlegroups.com.

[deal.II] Code tries to access DoF indices on artificial cells, even though those cells should not be artificial

['Maxi Miller' via deal.II User Group]

I am trying to port example 15 to a MPI-supported version, using the 
geometric multigrid method as preconditioner (to get more familiar with the 
implementation). I followed example 50 while doing that, but currently I 
get the error
An error occurred in line <3866> of file 
<~/Downloads/git-files/dealii/include/deal.II/dofs/dof_accessor.templates.h> 
in function
void dealii::DoFCellAccessor::get_dof_values(const 
InputVector&, ForwardIterator, ForwardIterator) const [with InputVector = 
dealii::TrilinosWrappers::MPI::Vector; ForwardIterator = double*; 
DoFHandlerType = dealii::DoFHandler<2, 2>; bool level_dof_access = true]
The violated condition was: 
this->is_artificial() == false
Additional information: 
Can't ask for DoF indices on artificial cells.
as soon as I try to run my program using several MPI threads. I tried to 
guard that part using if (cell->level_subdomain_id() == 
triangulation.locally_owned_subdomain()) as suggested in step-50, but that 
did not help. Why am I still trying to request data from artificial cells, 
even though I guarded that code part (at least I assume it guards it 
accordingly).
Thanks!

-- 
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/d7eecf00-5fce-46c7-a191-a6d6261dac35%40googlegroups.com.
/* -
 *
 * Copyright (C) 2012 - 2019 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, redistribute
 * it, and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * The full text of the license can be found in the file LICENSE.md at
 * the top level directory of deal.II.
 *
 * -

 *
 * Author: Sven Wetterauer, University of Heidelberg, 2012
 */


// @sect3{Include files}

// The first few files have already been covered in previous examples and will
// thus not be further commented on.
#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 
#include 

#include 
#include 

#include 
#include 
#include 
#include 

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

#include 
#include 
#include 

#include 

#include 
#include 

// We will use adaptive mesh refinement between Newton iterations. To do so,
// we need to be able to work with a solution on the new mesh, although it was
// computed on the old one. The SolutionTransfer class transfers the solution
// from the old to the new mesh:

#include 

// We then open a namespace for this program and import everything from the
// dealii namespace into it, as in previous programs:
namespace Step15
{
using namespace dealii;

namespace LA{
using namespace dealii::LinearAlgebraTrilinos;
}

template 
struct ScratchData
{
	ScratchData(const Mapping &  mapping,
const FiniteElement ,
const unsigned intquadrature_degree,
const UpdateFlags update_flags)
		: fe_values(mapping, fe, QGauss(quadrature_degree), update_flags)
	{}
	ScratchData(const ScratchData _data)
		: fe_values(scratch_data.fe_values.get_mapping(),
	scratch_data.fe_values.get_fe(),
	scratch_data.fe_values.get_quadrature(),
	scratch_data.fe_values.get_update_flags())
	{}
	FEValues fe_values;
};

struct CopyData
{
	unsigned int level;
	FullMatrix   cell_matrix;
	Vector   cell_rhs;
	std::vector local_dof_indices;
	template 
	void reinit(const Iterator , unsigned int dofs_per_cell)
	{
		cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
		cell_rhs.reinit(dofs_per_cell);
		local_dof_indices.resize(dofs_per_cell);
		cell->get_active_or_mg_dof_indices(local_dof_indices);
		level = cell->level();
	}
};

// @sect3{The MinimalSurfaceProblem class template}

// The class template is basically the same as in step-6.  Three additions
// are made:
// - There are two solution vectors, one for the Newton update
//   $\delta u^n$, and one for the current iterate $u^n$.
// - The setup_system function takes an argument that denotes
//   whether this is the first time it is called or not. The difference is
//   that the first time around we need to distribute the degrees of freedom
//   and set the solution vector for 

[deal.II] Create a LinearOperator with a function instead of a matrix

['Maxi Miller' via deal.II User Group]

I tried to follow the examples for setting up a LinearOperator (such as 
shown in Step-22 and Step-20), but instead of a matrix I wanted to provide 
my own function which should serve as vmult-function (i.e. a matrix-free 
linear operator). Is that possible? My initial test was the following code:

  class approximator : public Subscriptor 
{ 
public: 
   approximator(std::function&, 
Vector&)> 
residual_function, 
const Vector _point) : 
residual_function(residual_function)
{ 
u = evaluation_point; 
u.add(1); 
epsilon_val = u.l1_norm() * 1e-6; 
residual_function(u, residual_0); 
}; 

void vmult(Vector , const Vector ) const; 

private: 
std::function&, 
   Vector&)> residual_function; 
double epsilon_val; 
Vector u, residual_0, residual_1; 
};


const auto op_M = linear_operator(jacobian_approximation(std::bind(&
MinimalSurfaceProblem::compute_residual,
   this,
   std::
placeholders::_1,
   std::
placeholders::_2), 

 present_solution));

but that resulted in
error: use of deleted function ‘dealii::LinearOperator dealii::linear_operator(Matrix&&) [with Range = dealii::Vector
; Domain = dealii::Vector; Payload = dealii::internal::
LinearOperatorImplementation::EmptyPayload; Matrix = Step15::
jacobian_approximation;  = void]’
  389 | 
 present_solution));
Are there other ways?
Thanks!


-- 
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/1cf547d8-4d9b-4826-b0c8-85c2bd2d2f3f%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Implicit stepping methods from the TimeStepping-namespace for large sparse matrices

['Maxi Miller' via deal.II User Group]

If I understand it correctly, that is equivalent to solving Ax = b, with A 
and b given, and x returned?

Am Montag, 15. Juli 2019 22:03:46 UTC+2 schrieb Daniel Arndt:
>
> Maxi,
>
> The interface for TimeStepping::evolve_one_time_step is:
>
> virtual double
>   evolve_one_time_step(
>   std::vector>
>   & F,
>   std::vector   VectorType(const double, const double, const VectorType &)>> _inverse,
>   double t,
>   double delta_t,
>   VectorType & y);
>
> and does not depend on UMFPACK at all. You just need to provide a 
> function-type object that can be used for evaluating the inverse of the 
> Jacobians and the right-hand side.
> You might also want to have a look at the test base/time_stepping_01 (
> https://github.com/dealii/dealii/blob/master/tests/base/time_stepping_01.cc
> ).
>
> Best,
> Daniel
>
>
> Am Mo., 15. Juli 2019 um 15:04 Uhr schrieb 'Maxi Miller' via deal.II User 
> Group >:
>
>> As far as I understand all implicit time-stepping methods in the 
>> TimeStepping-namespace take a function which invert the jacobian matrix and 
>> the mass matrix (combined). This is done using a sparse solver (UMFPACK). 
>> Is it that still possible for a larger system (10kk DoFs), or do I have to 
>> resort to other methods, after inverting the sparse matrix will result in a 
>> dense matrix (which likely will overflow my memory)?
>>
>> -- 
>> 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/ccfe4a75-775d-471b-b69e-4a886c2e443d%40googlegroups.com
>>  
>> <https://groups.google.com/d/msgid/dealii/ccfe4a75-775d-471b-b69e-4a886c2e443d%40googlegroups.com?utm_medium=email_source=footer>
>> .
>> For more options, visit https://groups.google.com/d/optout.
>>
>

-- 
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/9937e20a-029f-4bc1-a030-31dbcff5c422%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Implicit stepping methods from the TimeStepping-namespace for large sparse matrices

['Maxi Miller' via deal.II User Group]

As far as I understand all implicit time-stepping methods in the 
TimeStepping-namespace take a function which invert the jacobian matrix and 
the mass matrix (combined). This is done using a sparse solver (UMFPACK). 
Is it that still possible for a larger system (10kk DoFs), or do I have to 
resort to other methods, after inverting the sparse matrix will result in a 
dense matrix (which likely will overflow my memory)?

-- 
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/ccfe4a75-775d-471b-b69e-4a886c2e443d%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Matrix-Vector-multiplication on deal.II-grid

['Maxi Miller' via deal.II User Group]

It is difficult to write it as a single integral. The operation is similar 
to the split-step fourier method, i.e. transforming the column vector f(r) 
once using
g(rho)=2\pi\int_0^\infty rf(r)J_0(2\pi\rho r)dr,
multiplying it with a vector, and transforming it back using
f(r) = 2\pi\int_0^\infty\rho g(\rho)J_0(2\pi\rho r)d\rho
The operation is for radially symmetric systems, i.e. with z along the 
x-axis, and r along the y-axis. When starting on the left border with f_0, 
i.e. at position z = 0, doing the operation mentioned above gives the 
values for the nodes at z = 1, when enumerating the nodes from 0 to n along 
the z axis, and having equidistant nodes along the z-axis. Those integrals 
can be replaced by a matrix-vector-multiplication, thus making it easier to 
implement numerically.
Does that make sense?

Am Freitag, 12. Juli 2019 18:24:22 UTC+2 schrieb Daniel Arndt:
>
> Maxi,
>
> can you clarify the operator evaluation you want to perform in 
> mathematical terms (maybe as an integral)?
>
> Best,
> Daniel
>

-- 
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/da921bbb-aad9-4ee9-af37-24971c981016%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Matrix-Vector-multiplication on deal.II-grid

['Maxi Miller' via deal.II User Group]

Is it possible to "mis-"use the facilities of deal.II for calculating a NLS 
using the discrete hankel transformation? The idea here was to take an 
evenly refined grid, and calculate the values of each node column based on 
the values of the node column to the left/right times a constant matrix, 
while beginning on the right/left side, and iterating to the other side? 
That would simplify the integration of my equations into my FEM-simulation 
significantly, especially because I can use deal.II-facilities for creating 
and managing the grid.
Thanks!

-- 
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/2b224489-2826-4061-871b-c0ae89fbd070%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Deal.II-Installation on Cluster fails

['Maxi Miller' via deal.II User Group]

Yes, that works, thanks!


Am Donnerstag, 11. Juli 2019 17:01:21 UTC+2 schrieb Daniel Arndt:
>
> Maxi,
>
> Can you just try again with a recent commit?
>
> Best,
> Daniel
>

-- 
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/47daa518-ed5d-4806-b62c-0e1654137c8e%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Elementwise multiplication of Tensors

['Maxi Miller' via deal.II User Group]

Something like this: https://github.com/dealii/dealii/pull/8307?

Am Donnerstag, 6. Juni 2019 00:38:59 UTC+2 schrieb Wolfgang Bangerth:
>
> On 6/5/19 3:22 PM, 'Maxi Miller' via deal.II User Group wrote: 
> > is there a function for the elementwise multiplication of Tensors? I 
> could 
> > only find functions for the scalar product, but not the schur product. 
>
> You mean to compute 
>C_ij = A_ij * B_ij 
> ? No, there isn't one right now. But it should easy to implement. In 
> analogy 
> to scalar_product(), the best approach would probably a global function 
> schur_product(). Want to give this a try? We'd of course take the patch! 
>
> Cheers 
>   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/985bb2a0-3408-4c6e-9e89-89ae8b7fd090%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Elementwise multiplication of Tensors

['Maxi Miller' via deal.II User Group]

Hei,
is there a function for the elementwise multiplication of Tensors? I could 
only find functions for the scalar product, but not the schur product.
Thanks!

-- 
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/734ec1b2-bee8-45bf-ac5b-2c8e3cc77cba%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Interpolation between FE_Q- and FE_Bernstein elements

['Maxi Miller' via deal.II User Group]

On a second test I found out that by resetting the out-vector to 0 I get 
the expected output, regardless of direction. Does that make sense?

Am Donnerstag, 30. Mai 2019 22:25:59 UTC+2 schrieb Maxi Miller:
>
> That means that I don't have to create that matrix myself, but rather use 
> the matrix I get from FE_Q_*::get_interpolation_matrix()? 
> The most interesting bug my code has at the moment is:
> Q(a)->B(b): I get identical vectors, i.e. a == b
> B(b)->Q(a): Vectors are not identical anymore, a != b
> Q(a)->B(b): Even though that worked in the first round, it does not in the 
> second, a != b
> I do not understand that behavior yet.
>
> Am Donnerstag, 30. Mai 2019 19:00:02 UTC+2 schrieb Wolfgang Bangerth:
>>
>> On 5/30/19 9:06 AM, 'Maxi Miller' via deal.II User Group wrote: 
>> > I wanted to write some functions to convert my solution from using 
>> > FE_Q-elements to FE_Bernstein-elements (and back). For that I wrote two 
>> > functions, convert_feQ_to_feB() and convert_feB_to_feQ() (as listed in 
>> the 
>> > attachment). When converting from FE_Q-elements to 
>> FE_Bernstein-elements, I 
>> > get the same values for input and output, but on the way back I get 
>> wrong 
>> > values (currently I am just using a vector filled with the value 1), 
>> but I can 
>> > not find the reason why the conversion back is failing. Did I forget 
>> something 
>> > in the corresponding function? 
>>
>> I have to admit that I don't have the time right now to look at your 
>> code, but 
>> it's really just an interpolation problem to go from FE_B to FE_Q: you 
>> need to 
>> evaluate each shape function of the FE_B at the interpolation (=support) 
>> points of the FE_Q. This gives you a cell-local matrix of the form 
>>B_{ij} = \varphi_{Bernstein,j)(x_{j,Q}) 
>> (or with indices reversed -- check with the documentation). 
>>
>> By the way, the place to implement this is in the 
>> FE_Q_Base::get_interpolation_matrix() function for B -> Q (see 
>> source/fe/fe_base.cc around line 500), and maybe 
>> FE_Q_Bernstein::get_interpolation_matrix() if you come up with a clever 
>> way to 
>> define the operation. I know that the latter function says that you can't 
>> do 
>> it, but I suspect that if you at least define it in a way so that the 
>> interpolation of a function onto itself is the identity operation, and if 
>> the 
>> interpolation of a FE_Q function onto FE_Bernstein yields the same, 
>> continuous 
>> function, then it might still go into that place. 
>>
>> 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/fe785f12-beee-4f10-827d-5869c34d09ff%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
#include 

#include 
#include 
#include 

#include 
#include 
#include 

#include 

#include 
#include 
#include 

#include 

#include 
#include 

using namespace dealii;

template 
void convert_feB_to_feQ(const DoFHandler ,
		const DoFHandler ,
		const Vector ,
		Vector ){
	out = 0.;
	AssertDimension(in.size(), dhB.n_dofs());
	const FiniteElement  = dhB.get_fe();
	const FiniteElement  = dhq.get_fe();

	const ComponentMask fe_mask(feB.get_nonzero_components(0).size(), true);

	AssertDimension(fe_mask.size(), feB.get_nonzero_components(0).size());

	std::vector fe_to_real(fe_mask.size(),
		 numbers::invalid_unsigned_int);
	unsigned int  size = 0;
	for (unsigned int i = 0; i < fe_mask.size(); ++i)
	{
		if (fe_mask[i])
			fe_to_real[i] = size++;
	}

	const ComponentMask maskq(dim, true);

	FullMatrix transfer(feB.dofs_per_cell, feq.dofs_per_cell);
	FullMatrix local_transfer(feq.dofs_per_cell);
	const std::vector>  = feq.get_unit_support_points();

	std::vector fe_to_feq(feB.dofs_per_cell,
		numbers::invalid_unsigned_int);
	unsigned int  index = 0;
	for (unsigned int i = 0; i < feB.dofs_per_cell; ++i)
		if (fe_mask[feB.system_to_component_index(i).first])
			fe_to_feq[i] = index++;

	Assert(index == feB.dofs_per

Re: [deal.II] Interpolation between FE_Q- and FE_Bernstein elements

['Maxi Miller' via deal.II User Group]

That means that I don't have to create that matrix myself, but rather use 
the matrix I get from FE_Q_*::get_interpolation_matrix()? 
The most interesting bug my code has at the moment is:
Q(a)->B(b): I get identical vectors, i.e. a == b
B(b)->Q(a): Vectors are not identical anymore, a != b
Q(a)->B(b): Even though that worked in the first round, it does not in the 
second, a != b
I do not understand that behavior yet.

Am Donnerstag, 30. Mai 2019 19:00:02 UTC+2 schrieb Wolfgang Bangerth:
>
> On 5/30/19 9:06 AM, 'Maxi Miller' via deal.II User Group wrote: 
> > I wanted to write some functions to convert my solution from using 
> > FE_Q-elements to FE_Bernstein-elements (and back). For that I wrote two 
> > functions, convert_feQ_to_feB() and convert_feB_to_feQ() (as listed in 
> the 
> > attachment). When converting from FE_Q-elements to 
> FE_Bernstein-elements, I 
> > get the same values for input and output, but on the way back I get 
> wrong 
> > values (currently I am just using a vector filled with the value 1), but 
> I can 
> > not find the reason why the conversion back is failing. Did I forget 
> something 
> > in the corresponding function? 
>
> I have to admit that I don't have the time right now to look at your code, 
> but 
> it's really just an interpolation problem to go from FE_B to FE_Q: you 
> need to 
> evaluate each shape function of the FE_B at the interpolation (=support) 
> points of the FE_Q. This gives you a cell-local matrix of the form 
>B_{ij} = \varphi_{Bernstein,j)(x_{j,Q}) 
> (or with indices reversed -- check with the documentation). 
>
> By the way, the place to implement this is in the 
> FE_Q_Base::get_interpolation_matrix() function for B -> Q (see 
> source/fe/fe_base.cc around line 500), and maybe 
> FE_Q_Bernstein::get_interpolation_matrix() if you come up with a clever 
> way to 
> define the operation. I know that the latter function says that you can't 
> do 
> it, but I suspect that if you at least define it in a way so that the 
> interpolation of a function onto itself is the identity operation, and if 
> the 
> interpolation of a FE_Q function onto FE_Bernstein yields the same, 
> continuous 
> function, then it might still go into that place. 
>
> 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/b7102277-a66f-4e38-84a8-cd90326dca0a%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Interpolation between FE_Q- and FE_Bernstein elements

['Maxi Miller' via deal.II User Group]

I wanted to write some functions to convert my solution from using 
FE_Q-elements to FE_Bernstein-elements (and back). For that I wrote two 
functions, convert_feQ_to_feB() and convert_feB_to_feQ() (as listed in the 
attachment). When converting from FE_Q-elements to FE_Bernstein-elements, I 
get the same values for input and output, but on the way back I get wrong 
values (currently I am just using a vector filled with the value 1), but I 
can not find the reason why the conversion back is failing. Did I forget 
something in the corresponding function?
Thanks!

-- 
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/e288a022-ff5b-4876-8ca3-9807605ef2a3%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
#include 

#include 
#include 
#include 

#include 
#include 
#include 

#include 

#include 
#include 
#include 

#include 

#include 
#include 

using namespace dealii;

template 
void convert_feB_to_feQ(const DoFHandler ,
		const DoFHandler ,
		const Vector ,
		Vector ){
	AssertDimension(in.size(), dhB.n_dofs());
	const FiniteElement  = dhB.get_fe();
	const FiniteElement  = dhq.get_fe();

	const ComponentMask fe_mask(feB.get_nonzero_components(0).size(), true);

	AssertDimension(fe_mask.size(), feB.get_nonzero_components(0).size());

	std::vector fe_to_real(fe_mask.size(),
		 numbers::invalid_unsigned_int);
	unsigned int  size = 0;
	for (unsigned int i = 0; i < fe_mask.size(); ++i)
	{
		if (fe_mask[i])
			fe_to_real[i] = size++;
	}

	const ComponentMask maskq(dim, true);

	FullMatrix transfer(feq.dofs_per_cell, feB.dofs_per_cell);
	FullMatrix local_transfer(feB.dofs_per_cell);
	const std::vector>  = feq.get_unit_support_points();

	std::vector fe_to_feq(feB.dofs_per_cell,
		numbers::invalid_unsigned_int);
	unsigned int  index = 0;
	for (unsigned int i = 0; i < feB.dofs_per_cell; ++i)
		if (fe_mask[feB.system_to_component_index(i).first])
			fe_to_feq[i] = index++;

	Assert(index == feB.dofs_per_cell, ExcNotImplemented());

	for (unsigned int j = 0; j < feB.dofs_per_cell; ++j)
	{
		const unsigned int comp_j = feB.system_to_component_index(j).first;
		if (fe_mask[comp_j])
			for (unsigned int i = 0; i < points.size(); ++i)
			{
if (fe_to_real[comp_j] ==
		feB.system_to_component_index(i).first)
	local_transfer(i, fe_to_feq[j]) =
			feq.shape_value(j, points[i]);
			}
	}

	local_transfer.invert(local_transfer);
	for (unsigned int i = 0; i < feq.dofs_per_cell; ++i)
		if (fe_to_feq[i] != numbers::invalid_unsigned_int)
			for (unsigned int j = 0; j < feB.dofs_per_cell; ++j)
transfer(i, j) = local_transfer(fe_to_feq[i], j);

	VectorTools::interpolate(dhB, dhq, transfer, in, out);
}

template 
void convert_feQ_to_feB(const DoFHandler ,
		const DoFHandler ,
		const Vector ,
		Vector ){

	AssertDimension(in.size(), dhq.n_dofs());
	const FiniteElement  = dhq.get_fe();
	const FiniteElement  = dhb.get_fe();

	const ComponentMask fe_mask(feq.get_nonzero_components(0).size(), true);

	AssertDimension(fe_mask.size(), feq.get_nonzero_components(0).size());

	std::vector fe_to_real(fe_mask.size(),
		 numbers::invalid_unsigned_int);
	unsigned int  size = 0;
	for (unsigned int i = 0; i < fe_mask.size(); ++i)
	{
		if (fe_mask[i])
			fe_to_real[i] = size++;
	}

	const ComponentMask maskq(dim, true);

	FullMatrix transfer(feq.dofs_per_cell, feB.dofs_per_cell);
	FullMatrix local_transfer(feB.dofs_per_cell);
	const std::vector>  = feq.get_unit_support_points();


	std::vector fe_to_feq(feq.dofs_per_cell,
		numbers::invalid_unsigned_int);
	unsigned int  index = 0;
	for (unsigned int i = 0; i < feq.dofs_per_cell; ++i)
		if (fe_mask[feq.system_to_component_index(i).first])
			fe_to_feq[i] = index++;

	Assert(index == feq.dofs_per_cell, ExcNotImplemented());

	for (unsigned int j = 0; j < feq.dofs_per_cell; ++j)
	{
		const unsigned int comp_j = feq.system_to_component_index(j).first;
		if (fe_mask[comp_j])
			for (unsigned int i = 0; i < points.size(); ++i)
			{
if (fe_to_real[comp_j] ==
		feq.system_to_component_index(i).first)
	local_transfer(i, fe_to_feq[j]) =
			feq.shape_value(j, points[i]);
			}
	}

	local_transfer.invert(local_transfer);
	for (unsigned int i = 0; i < feq.dofs_per_cell; ++i)
		if (fe_to_feq[i] != numbers::invalid_unsigned_int)
			for (unsigned int j = 0; j < feB.dofs_per_cell; ++j)
transfer(i, j) = local_transfer(fe_to_feq[i], j);

	VectorTools::interpolate(dhq, dhb, transfer, in, out);
}

template 
void test_interpolation(const FiniteElement _1,

Re: [deal.II] prepare_for_coarsening_and_refinement() for FE_Bernstein-elements

['Maxi Miller' via deal.II User Group]

ry, but are there any obvious bugs in 
the code?
Thanks!

Am Dienstag, 21. Mai 2019 14:08:32 UTC+2 schrieb Maxi Miller:
>
> According to 
> https://www.dealii.org/developer/doxygen/deal.II/classFE__Bernstein.html#ad881d4e04f5e699e4a8bda5b4e9f5265
>  
> that is impossible. Would an alternative solution be to interpolate 
> everything onto a grid using FE_Q-elements, transfer the result, and 
> interpolate the results back?
>
> Am Dienstag, 21. Mai 2019 03:41:58 UTC+2 schrieb Wolfgang Bangerth:
>>
>> On 5/20/19 5:18 AM, 'Maxi Miller' via deal.II User Group wrote: 
>> > Trying to transfer solutions from an old grid to a new grid with 
>> > parallel::distributed::SolutionTransfer> TrilinosWrappers::MPI::Vector, 
>> > DoFHandler> solution_transfer() and 
>> > solution_transfer.prepare_for_coarsening_and_refinement() fails with 
>> > | 
>> > Anerror occurred inline <71>of file 
>> > <~/Downloads/git-files/dealii/source/fe/fe_bernstein.cc>infunction 
>> > 
>> constdealii::FullMatrix::FE_Bernstein::get_restriction_matrix(unsignedint,constdealii::RefinementCase&)const[withintdim
>>  
>>
>> > =2;intspacedim =2] 
>> > Theviolated condition was: 
>> > false 
>> > Additionalinformation: 
>> > Youare trying to access the matrices that describe how to restrict a 
>> finite 
>> > element functionfromthe children of one cell to the finite element 
>> space 
>> > definedon their parent (i.e.,the 
>> > 'restriction'or'projection'matrices).However,the current finite element 
>> can 
>> > either notdefine thissort of operation,orit has notyet been 
>> implemented. 
>> > | 
>> > 
>> > 
>> > Which function do I have to modify to make that work? 
>>
>> You will have to implement the restriction matrices for the Bernstein 
>> element. 
>> These are usually set in the constructor of the finite element classes. 
>> Take a 
>> look at the other FE classes for examples. 
>>
>> The error you see is simply because you are accessing a field that has 
>> never 
>> been initialized. 
>>
>> 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/fa00b08a-8f8f-44d4-bfba-16d0c3735cf4%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] prepare_for_coarsening_and_refinement() for FE_Bernstein-elements

['Maxi Miller' via deal.II User Group]

According to 
https://www.dealii.org/developer/doxygen/deal.II/classFE__Bernstein.html#ad881d4e04f5e699e4a8bda5b4e9f5265
 
that is impossible. Would an alternative solution be to interpolate 
everything onto a grid using FE_Q-elements, transfer the result, and 
interpolate the results back?

Am Dienstag, 21. Mai 2019 03:41:58 UTC+2 schrieb Wolfgang Bangerth:
>
> On 5/20/19 5:18 AM, 'Maxi Miller' via deal.II User Group wrote: 
> > Trying to transfer solutions from an old grid to a new grid with 
> > parallel::distributed::SolutionTransfer TrilinosWrappers::MPI::Vector, 
> > DoFHandler> solution_transfer() and 
> > solution_transfer.prepare_for_coarsening_and_refinement() fails with 
> > | 
> > Anerror occurred inline <71>of file 
> > <~/Downloads/git-files/dealii/source/fe/fe_bernstein.cc>infunction 
> > 
> constdealii::FullMatrix::FE_Bernstein::get_restriction_matrix(unsignedint,constdealii::RefinementCase&)const[withintdim
>  
>
> > =2;intspacedim =2] 
> > Theviolated condition was: 
> > false 
> > Additionalinformation: 
> > Youare trying to access the matrices that describe how to restrict a 
> finite 
> > element functionfromthe children of one cell to the finite element space 
> > definedon their parent (i.e.,the 
> > 'restriction'or'projection'matrices).However,the current finite element 
> can 
> > either notdefine thissort of operation,orit has notyet been implemented. 
> > | 
> > 
> > 
> > Which function do I have to modify to make that work? 
>
> You will have to implement the restriction matrices for the Bernstein 
> element. 
> These are usually set in the constructor of the finite element classes. 
> Take a 
> look at the other FE classes for examples. 
>
> The error you see is simply because you are accessing a field that has 
> never 
> been initialized. 
>
> 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/cdf53472-6b18-4c33-b5a4-513c4025da08%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Interpolate correct solution onto mesh using interpolate() with FE_Bernstein elements

['Maxi Miller' via deal.II User Group]

I wrote a drop-in function which can replace VectorTools::interpolate, but 
works for both FE_Q and FE_Bernstein (not tested yet for other elements). 
Is it useful to have something like that in the library, and if yes, what 
is the best way for proof-checking and integrating?
Thanks!

Am Donnerstag, 16. Mai 2019 22:00:23 UTC+2 schrieb luca.heltai:
>
>
> > On 16 May 2019, at 9:16, 'Maxi Miller' via deal.II User Group <
> dea...@googlegroups.com > wrote: 
> > 
> > I would like to use the interpolate()-function to interpolate the 
> correct function onto my mesh for comparison during some simple tests. This 
> works fine when using FE_Q-elements, but fails when using FE_Bernstein 
> elements. 
>
> FE_Bernstein, like BSplines, are not interpolatory, there is not a single 
> definition of “interpolation”, since there are no support points. 
>
> If you want your function to be interpolated at some points, then you’d 
> have to construct the interpolation matrix yourself, by solving in some 
> sense the following system 
>
> sum_i  B_i(x_j) u_i = f(x_j) 
>
> provided that you specify exactly the same number of interpolation points 
> as there are degrees of freedom. 
>
> One easy way to do this, is to interpolate to a dofhandler with FE_Q(d), 
> with the same degree of the FEBernstein, and then call 
>
>
> void VectorTools::interpolate(const DoFHandler< dim, 
> spacedim > & dof_1, 
> const DoFHandler< dim, spacedim > & dof_2, 
> const FullMatrix< double > & transfer, 
> const InVector & data_1, 
> OutVector & data_2 
> ) 
>
> where the FullMatrix< double > & transfer represents the local 
> interpolation between FE_Bernstein and FE_Q. 
>
>
> https://www.dealii.org/current/doxygen/deal.II/namespaceVectorTools.html#a5e3af70a47cedfaf361cf5c621e94e3d
>  
>
> This is used, for example, here: 
>
>
> https://github.com/dealii/dealii/blob/master/include/deal.II/numerics/vector_tools.templates.h#L9120
>  
>
> You can do the exact same thing. Construct a local interpolation matrix, 
> and then call the interpolate function above passing this matrix (and the 
> vector you interpolated on the standard FE_Q space). 
>
>
> https://github.com/dealii/dealii/blob/master/include/deal.II/numerics/vector_tools.templates.h#L9208
>  
>
> Best, 
> Luca. 
>
>
>
> > The error message I get is 
> > An error occurred in line <3098> of file 
> <~/Downloads/git-files/dealii/source/fe/fe_values.cc> in function 
> > dealii::FEValuesBase::FEValuesBase(unsigned int, 
> unsigned int, dealii::UpdateFlags, const dealii::Mapping&, 
> const dealii::FiniteElement&) [with int dim = 2; int 
> spacedim = 2] 
> > The violated condition was: 
> > n_q_points > 0 
> > Additional information: 
> > There is nothing useful you can do with an FEValues object when 
> using a quadrature formula with zero quadrature points! 
> > 
> > 
> > Is there an alternative way to interpolate a function onto the mesh when 
> using FE_Bernstein elements? 
> > 
> > -- 
> > 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/2f5b238a-b78a-4097-b492-2e0ebd82ec93%40googlegroups.com.
>  
>
> > For more options, visit https://groups.google.com/d/optout. 
>
>

-- 
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/f8829f33-a47a-41fd-a74e-159b8ac698c0%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Interpolate correct solution onto mesh using interpolate() with FE_Bernstein elements

['Maxi Miller' via deal.II User Group]

That works, thanks!

Am Donnerstag, 16. Mai 2019 22:00:23 UTC+2 schrieb luca.heltai:
>
>
> > On 16 May 2019, at 9:16, 'Maxi Miller' via deal.II User Group <
> dea...@googlegroups.com > wrote: 
> > 
> > I would like to use the interpolate()-function to interpolate the 
> correct function onto my mesh for comparison during some simple tests. This 
> works fine when using FE_Q-elements, but fails when using FE_Bernstein 
> elements. 
>
> FE_Bernstein, like BSplines, are not interpolatory, there is not a single 
> definition of “interpolation”, since there are no support points. 
>
> If you want your function to be interpolated at some points, then you’d 
> have to construct the interpolation matrix yourself, by solving in some 
> sense the following system 
>
> sum_i  B_i(x_j) u_i = f(x_j) 
>
> provided that you specify exactly the same number of interpolation points 
> as there are degrees of freedom. 
>
> One easy way to do this, is to interpolate to a dofhandler with FE_Q(d), 
> with the same degree of the FEBernstein, and then call 
>
>
> void VectorTools::interpolate(const DoFHandler< dim, 
> spacedim > & dof_1, 
> const DoFHandler< dim, spacedim > & dof_2, 
> const FullMatrix< double > & transfer, 
> const InVector & data_1, 
> OutVector & data_2 
> ) 
>
> where the FullMatrix< double > & transfer represents the local 
> interpolation between FE_Bernstein and FE_Q. 
>
>
> https://www.dealii.org/current/doxygen/deal.II/namespaceVectorTools.html#a5e3af70a47cedfaf361cf5c621e94e3d
>  
>
> This is used, for example, here: 
>
>
> https://github.com/dealii/dealii/blob/master/include/deal.II/numerics/vector_tools.templates.h#L9120
>  
>
> You can do the exact same thing. Construct a local interpolation matrix, 
> and then call the interpolate function above passing this matrix (and the 
> vector you interpolated on the standard FE_Q space). 
>
>
> https://github.com/dealii/dealii/blob/master/include/deal.II/numerics/vector_tools.templates.h#L9208
>  
>
> Best, 
> Luca. 
>
>
>
> > The error message I get is 
> > An error occurred in line <3098> of file 
> <~/Downloads/git-files/dealii/source/fe/fe_values.cc> in function 
> > dealii::FEValuesBase::FEValuesBase(unsigned int, 
> unsigned int, dealii::UpdateFlags, const dealii::Mapping&, 
> const dealii::FiniteElement&) [with int dim = 2; int 
> spacedim = 2] 
> > The violated condition was: 
> > n_q_points > 0 
> > Additional information: 
> > There is nothing useful you can do with an FEValues object when 
> using a quadrature formula with zero quadrature points! 
> > 
> > 
> > Is there an alternative way to interpolate a function onto the mesh when 
> using FE_Bernstein elements? 
> > 
> > -- 
> > 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/2f5b238a-b78a-4097-b492-2e0ebd82ec93%40googlegroups.com.
>  
>
> > For more options, visit https://groups.google.com/d/optout. 
>
>

-- 
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/054854ed-da4c-4e54-8270-bf8b658e9798%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] prepare_for_coarsening_and_refinement() for FE_Bernstein-elements

['Maxi Miller' via deal.II User Group]

Trying to transfer solutions from an old grid to a new grid with 
parallel::distributed::SolutionTransfer> solution_transfer() and 
solution_transfer.prepare_for_coarsening_and_refinement() fails with 
An error occurred in line <71> of file <~/Downloads/git-files/dealii/source/
fe/fe_bernstein.cc> in function
const dealii::FullMatrix& dealii::FE_Bernstein::
get_restriction_matrix(unsigned int, const dealii::RefinementCase&) 
const [with int dim = 2; int spacedim = 2]
The violated condition was: 
false
Additional information: 
You are trying to access the matrices that describe how to restrict a 
finite element function from the children of one cell to the finite element 
space defined on their parent (i.e., the 'restriction' or 'projection' 
matrices). However, the current finite element can either not define this 
sort of operation, or it has not yet been implemented.


Which function do I have to modify to make that work?
Thanks!

-- 
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/1489fb1a-25f1-4eaa-9cb4-69b4964cf7ba%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Applying boundary masks when using project_boundary_values()

['Maxi Miller' via deal.II User Group]

If I understand that correctly, the correct order for doing that would be:
-> Interpolate boundary condition without any component mask, giving me a 
vector
-> Loop over all elements in the vector, and cross-check with the mask 
(using extract_dofs()) if they should be interpolated in that way
-> If yes, copy them to a second vector
-> Repeat with the next boundary condition and new component mask
-> After finishing: Replace first vector with second vector (or use a 
temporary vector as first vector)
Is that correct?
In addition: How do I treat different boundary markings, i.e. 
dirichlet-conditions on 0, neumann-conditions on 1, and similar stuff?
Thanks!

Am Donnerstag, 16. Mai 2019 19:20:59 UTC+2 schrieb Wolfgang Bangerth:
>
> On 5/15/19 8:37 AM, 'Maxi Miller' via deal.II User Group wrote: 
> > I have to replace VectorTools::interpolate_boundary_values with 
> > VectorTools::project_boundary_values at several places in my project. I 
> could 
> > apply the boundary values only for certain parts in my FESystem using 
> > component masks, but I can not use those masks in 
> project_boundary_values(). 
> > What would be the best approach for replacing them in this function? 
>
> That's a question that has been asked repeatedly in the recent past. In 
> other 
> words, there is clearly a demand to make this work, but nobody has found 
> the 
> time to implement it. So, if you want to do it right, the best place to 
> implement this is probably in the library itself. 
>
> But if you want to do it cheaply, you can always interpolate *all* 
> components 
> of your solution. You'd then either get a std::map or AffineMatrix object, 
> and 
> you will have to filter out those elements that do not match the 
> components 
> you are interested in. In practice, it's probably easier to just go 
> through 
> all elements of these objects, see which ones you care about (using 
> DoFTools::extracy_dofs()), and copy them to a second object. 
>
> If you want help to get this into the library, let us know! 
>
> 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/0ed11fe9-4ee1-4b91-a86f-5db7b696c419%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Interpolate correct solution onto mesh using interpolate() with FE_Bernstein elements

['Maxi Miller' via deal.II User Group]

I will take a look, thanks for the help! Would it be useful to have it in a 
form which could be introduced in the library later (as drop-in replacement 
for interpolate(), for example)?

Am Donnerstag, 16. Mai 2019 22:00:23 UTC+2 schrieb luca.heltai:
>
>
> > On 16 May 2019, at 9:16, 'Maxi Miller' via deal.II User Group <
> dea...@googlegroups.com > wrote: 
> > 
> > I would like to use the interpolate()-function to interpolate the 
> correct function onto my mesh for comparison during some simple tests. This 
> works fine when using FE_Q-elements, but fails when using FE_Bernstein 
> elements. 
>
> FE_Bernstein, like BSplines, are not interpolatory, there is not a single 
> definition of “interpolation”, since there are no support points. 
>
> If you want your function to be interpolated at some points, then you’d 
> have to construct the interpolation matrix yourself, by solving in some 
> sense the following system 
>
> sum_i  B_i(x_j) u_i = f(x_j) 
>
> provided that you specify exactly the same number of interpolation points 
> as there are degrees of freedom. 
>
> One easy way to do this, is to interpolate to a dofhandler with FE_Q(d), 
> with the same degree of the FEBernstein, and then call 
>
>
> void VectorTools::interpolate(const DoFHandler< dim, 
> spacedim > & dof_1, 
> const DoFHandler< dim, spacedim > & dof_2, 
> const FullMatrix< double > & transfer, 
> const InVector & data_1, 
> OutVector & data_2 
> ) 
>
> where the FullMatrix< double > & transfer represents the local 
> interpolation between FE_Bernstein and FE_Q. 
>
>
> https://www.dealii.org/current/doxygen/deal.II/namespaceVectorTools.html#a5e3af70a47cedfaf361cf5c621e94e3d
>  
>
> This is used, for example, here: 
>
>
> https://github.com/dealii/dealii/blob/master/include/deal.II/numerics/vector_tools.templates.h#L9120
>  
>
> You can do the exact same thing. Construct a local interpolation matrix, 
> and then call the interpolate function above passing this matrix (and the 
> vector you interpolated on the standard FE_Q space). 
>
>
> https://github.com/dealii/dealii/blob/master/include/deal.II/numerics/vector_tools.templates.h#L9208
>  
>
> Best, 
> Luca. 
>
>
>
> > The error message I get is 
> > An error occurred in line <3098> of file 
> <~/Downloads/git-files/dealii/source/fe/fe_values.cc> in function 
> > dealii::FEValuesBase::FEValuesBase(unsigned int, 
> unsigned int, dealii::UpdateFlags, const dealii::Mapping&, 
> const dealii::FiniteElement&) [with int dim = 2; int 
> spacedim = 2] 
> > The violated condition was: 
> > n_q_points > 0 
> > Additional information: 
> > There is nothing useful you can do with an FEValues object when 
> using a quadrature formula with zero quadrature points! 
> > 
> > 
> > Is there an alternative way to interpolate a function onto the mesh when 
> using FE_Bernstein elements? 
> > 
> > -- 
> > 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/2f5b238a-b78a-4097-b492-2e0ebd82ec93%40googlegroups.com.
>  
>
> > For more options, visit https://groups.google.com/d/optout. 
>
>

-- 
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/592f1e2a-40db-47ea-885d-e3b5a1701bd1%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Interpolate correct solution onto mesh using interpolate() with FE_Bernstein elements

['Maxi Miller' via deal.II User Group]

I would like to use the interpolate()-function to interpolate the correct 
function onto my mesh for comparison during some simple tests. This works 
fine when using FE_Q-elements, but fails when using FE_Bernstein elements. 
The error message I get is 
An error occurred in line <3098> of file <~/Downloads/git-files/dealii/
source/fe/fe_values.cc> in function
dealii::FEValuesBase::FEValuesBase(unsigned int, unsigned 
int, dealii::UpdateFlags, const dealii::Mapping&, const 
dealii::FiniteElement&) [with int dim = 2; int spacedim = 2]
The violated condition was: 
n_q_points > 0
Additional information: 
There is nothing useful you can do with an FEValues object when using a 
quadrature formula with zero quadrature points!


Is there an alternative way to interpolate a function onto the mesh when 
using FE_Bernstein elements?

-- 
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/2f5b238a-b78a-4097-b492-2e0ebd82ec93%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Applying boundary masks when using project_boundary_values()

['Maxi Miller' via deal.II User Group]

I have to replace VectorTools::interpolate_boundary_values with 
VectorTools::project_boundary_values at several places in my project. I 
could apply the boundary values only for certain parts in my FESystem using 
component masks, but I can not use those masks in 
project_boundary_values(). What would be the best approach for replacing 
them in this function?

-- 
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/a86d0863-20f1-43ca-8742-cefbd81a6db7%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Spline basìs function instead of gaussian basis function

['Maxi Miller' via deal.II User Group]

When using FE_Q-elements as basis function, Lagrange-polynoms are used, as 
far as I understand. Is it possible to replace those with B-splines as 
basis functions?

-- 
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/4786799e-c8ec-45c4-955b-c25f3580c066%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] scratch_data_03.cc from test/meshworker produces strange results when refining adaptively

['Maxi Miller' via deal.II User Group]

I wanted to test the new ScratchData- and CopyData-Classes for DG-subfaces, 
and thus used Test 3 with a locally refined grid. The result shows 
significant distortions at the transitions between refinement levels. 
Should that be, or is that a bug? Or is it just the missing handling of 
subfaces in that test? The demo-file is attached, including modifications

-- 
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.
/* -
 *
 * Copyright (C) 2018 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, redistribute
 * it, and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * The full text of the license can be found in the file LICENSE.md at
 * the top level directory of deal.II.
 *
 * -
 */

// Solve Laplacian using SIPG + mesh_loop + ScratchData

#include 
#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;
using namespace MeshWorker;

template 
void
test()
{
Triangulation tria;
FE_DGQfe(1);
DoFHandlerdh(tria);

FunctionParser rhs_function("1");
FunctionParser boundary_function("0");

AffineConstraints constraints;
constraints.close();

GridGenerator::hyper_cube(tria);
const size_t refinement_val = 5;
tria.refine_global(refinement_val);

tria.execute_coarsening_and_refinement();
#ifdef REFINE_ADAPTIVELY
for (unsigned int i = 0; i < 2; i++)
{
typename Triangulation::active_cell_iterator
cell = tria.begin_active(),
endc = tria.end();
for (; cell != endc; cell++)
{
if ((cell->center()[1]) > 0.75 )
{
if ((cell->center()[0] > 0.8)  || (cell->center()[0] < 0.2))
cell->set_refine_flag();
}
}
tria.execute_coarsening_and_refinement();
}
#endif

dh.distribute_dofs(fe);


SparsityPattern sparsity;

{
DynamicSparsityPattern dsp(dh.n_dofs(), dh.n_dofs());
DoFTools::make_flux_sparsity_pattern(dh, dsp);
sparsity.copy_from(dsp);
}

SparseMatrix matrix;
matrix.reinit(sparsity);

Vector solution(dh.n_dofs());
Vector rhs(dh.n_dofs());

QGauss quad(3);
QGauss face_quad(3);

UpdateFlags cell_flags = update_values | update_gradients |
update_quadrature_points | update_JxW_values;
UpdateFlags face_flags = update_values | update_gradients |
update_quadrature_points |
update_normal_vectors | update_JxW_values;

// Stabilization for SIPG
double gamma = 100;

using ScratchData = MeshWorker::ScratchData;
using CopyData = MeshWorker::CopyData<1 + GeometryInfo::faces_per_cell,
1,
1 + GeometryInfo::faces_per_cell>;

ScratchData scratch(fe, quad, cell_flags, face_quad, face_flags);
CopyDatacopy(fe.dofs_per_cell);

auto cell = dh.begin_active();
auto endc = dh.end();

typedef decltype(cell) Iterator;

auto cell_worker =
[_function](const Iterator , ScratchData , CopyData ) {
const auto  = s.reinit(cell);
const auto  = s.get_JxW_values();
const auto= s.get_quadrature_points();

c = 0;

c.local_dof_indices[0] = s.get_local_dof_indices();

for (unsigned int q = 0; q < p.size(); ++q)
for (unsigned int i = 0; i < fev.dofs_per_cell; ++i)
{
for (unsigned int j = 0; j < fev.dofs_per_cell; ++j)
{
c.matrices[0](i, j) +=
fev.shape_grad(i, q) * fev.shape_grad(j, q) * JxW[q];
}
c.vectors[0](i) +=
fev.shape_value(i, q) * rhs_function.value(p[q]) * JxW[q];
}
  };

auto boundary_worker = [gamma, _function](const Iterator &cell,
const unsigned int ,
ScratchData &   s,
CopyData &  c) {
const auto  = s.reinit(cell, f);
const auto  = s.get_JxW_values();
const auto= s.get_quadrature_points();
const auto= s.get_normal_vectors();

for (unsigned int q = 0; q < p.size(); ++q)
for (unsigned int i = 0; i < fev.dofs_per_cell; ++i)
{

Re: [deal.II] Re: Integrating a deal.II-specific function in a NOX-Interface for MPI-threads

['Maxi Miller' via deal.II User Group]

Am Dienstag, 30. April 2019 05:59:05 UTC+2 schrieb Wolfgang Bangerth:
>
>
> Maxi, 
>
> > What did you have to do? I was going to reproduce your problem today 
> by 
> > installing a Trilinos version that has NOX enabled. Is this moot 
> now, 
> > i.e., was it a bug in your code or did you just work around the 
> issue in 
> > some way that doesn't expose it? 
> > 
> > NOX expects vectors containing only the local elements, but no ghost 
> elements. 
> > Thus I had to initialize all vectors going in or out from any 
> NOX-related 
> > function using locally_owned_dofs, and copy accordingly if external 
> vectors 
> > contained ghost elements. 
>
> I see. So this is a NOX requirement, not something that we could have done 
> anything about on the deal.II side? 
>
> No, as far as I understand. NOX needs those vectors to be in the correct 
way, else it will not work. 

>
> > The solver does not converge, and the output looks as if it is using 
> > Dirichlet-conditions with u = 0, independently of the "real" boundary 
> conditions. 
>
> I don't know NOX, but is it using an update that it adds to the solution 
> in 
> each step? If so, you need to have the correct boundary conditions in the 
> initial guess already. What happens if you only allow NOX to do zero or 
> one 
> iterations? 
>
> Best 
>   W. 
>
> Zero iterations is not allowed, the first iteration already goes way off 
(expected highest value: 1, calculated value: 1e6), regardless if I set 
boundary conditions already to the solution vector I feed to NOX, or if I 
set the boundary conditions during the assembly of the right hand side. 

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

Re: [deal.II] Re: Integrating a deal.II-specific function in a NOX-Interface for MPI-threads

['Maxi Miller' via deal.II User Group]

Am Montag, 29. April 2019 22:36:55 UTC+2 schrieb Wolfgang Bangerth:
>
> On 4/29/19 7:35 AM, 'Maxi Miller' via deal.II User Group wrote: 
> > Could get it to work, now the code works equally fine for one or several 
> > MPI threads. 
>
> What did you have to do? I was going to reproduce your problem today by 
> installing a Trilinos version that has NOX enabled. Is this moot now, 
> i.e., was it a bug in your code or did you just work around the issue in 
> some way that doesn't expose it? 
>
> NOX expects vectors containing only the local elements, but no ghost 
elements. Thus I had to initialize all vectors going in or out from any 
NOX-related function using locally_owned_dofs, and copy accordingly if 
external vectors contained ghost elements. 

>
> > Still, the next open question is how to apply non-zero 
> > Dirichlet boundary conditions (as in the file, for cos(pi*x)*cos(pi*y). 
> > When applying them as I do in the file, I get the warning that the 
> > solver is not converged. When using Neumann-Conditions, the solver 
> > converges and I get the correct output. What would be the solution for 
> that? 
>
> The question to ask is who is right. When you get the message that the 
> solver is not converged, is this correct? That is, does it converge if 
> you allow more iterations? What if you use a case where you know the 
> exact solution -- is the computed solution actively wrong? 
>
The solver does not converge, and the output looks as if it is using 
Dirichlet-conditions with u = 0, independently of the "real" boundary 
conditions.  

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

Re: [deal.II] Re: Integrating a deal.II-specific function in a NOX-Interface for MPI-threads

['Maxi Miller' via deal.II User Group]

Could get it to work, now the code works equally fine for one or several 
MPI threads. Still, the next open question is how to apply non-zero 
Dirichlet boundary conditions (as in the file, for cos(pi*x)*cos(pi*y). 
When applying them as I do in the file, I get the warning that the solver 
is not converged. When using Neumann-Conditions, the solver converges and I 
get the correct output. What would be the solution for that?
Thanks!

Am Montag, 29. April 2019 10:20:13 UTC+2 schrieb Maxi Miller:
>
> Now with file...
>
> Am Montag, 29. April 2019 10:19:44 UTC+2 schrieb Maxi Miller:
>>
>> The functions computePreconditioner and computeJacobian are necessary, 
>> else I will get compilation problems. Thus I think the current version is 
>> the best MWE I can get. Running with one thread gives a "Test passed", 
>> using more threads will result in "Convergence failed".
>>
>> Am Samstag, 27. April 2019 05:43:47 UTC+2 schrieb Wolfgang Bangerth:
>>>
>>> On 4/26/19 2:28 PM, 'Maxi Miller' via deal.II User Group wrote: 
>>> > I think the current version is the smallest version I can get. It will 
>>> work 
>>> > for a single thread, but not for two or more threads. 
>>>
>>> I tried to run this on my machine, but I don't have NOX as part of my 
>>> Trilinos 
>>> installation. That might have to wait for next week. 
>>>
>>> But I think there is still room for making the program shorter. The 
>>> point is 
>>> that the parts of the program that run before the error happens do not 
>>> actually have to make sense. For example, could you replace your own 
>>> boundary 
>>> values class by ZeroFunction? Instead of assembling something real, 
>>> could you 
>>> just use the identity matrix on each cell? 
>>>
>>> Other things you don't need: Timers, convergence tables, etc. Which of 
>>> the 
>>> parameters you set in solve_for_NOX() are necessary? (Necessary to 
>>> reproduce 
>>> the bug; I know they're necessary for your application, but that's 
>>> irrelevant 
>>> here.) Could you replace building the residual by just putting random 
>>> stuff in 
>>> the vector? 
>>>
>>> Similarly, the computeJacobian and following functions really just check 
>>> conditions, but when you run the program to illustrate the bug you're 
>>> chasing, 
>>> do these checks trigger? If not, remove them, then remove the calls to 
>>> these 
>>> functions (because they don't do anything any more), then remove the 
>>> whole 
>>> function. 
>>>
>>> I think there's still room to make this program smaller by at least a 
>>> factor 
>>> of 2 or 3 or 4 :-) 
>>>
>>> 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.
For more options, visit https://groups.google.com/d/optout.
//Nox include files
#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 
#include 
#include 

#include 
#include 
#include 

#include 
#include 

#include 
#include 
#include 

#include 

#include 
#include 

#include 
#include 

// NOX Objects
#include "NOX.H"
#include "NOX_Epetra.H"

// Trilinos Objects
#ifdef HAVE_MPI
#include "Epetra_MpiComm.h"
#else
#include "Epetra_SerialComm.h"
#endif
#include "Epetra_Vector.h"
#include "Epetra_Operator.h"
#include "Epetra_RowMatrix.h"
#include "NOX_Epetra_Interface_Required.H" // base class
#include "NOX_Epetra_Interface_Jacobian.H" // base class
#include "NOX_Epetra_Interface_Preconditioner.H" // base class
#include "Epetra_CrsMatrix.h"
#include "Epetra_Map.h"
#include "Epetra_LinearProblem.h"
#include "AztecOO.h"
#include "Teuchos_StandardCatchMacros.hpp"

/

Re: [deal.II] Re: Integrating a deal.II-specific function in a NOX-Interface for MPI-threads

['Maxi Miller' via deal.II User Group]

Now with file...

Am Montag, 29. April 2019 10:19:44 UTC+2 schrieb Maxi Miller:
>
> The functions computePreconditioner and computeJacobian are necessary, 
> else I will get compilation problems. Thus I think the current version is 
> the best MWE I can get. Running with one thread gives a "Test passed", 
> using more threads will result in "Convergence failed".
>
> Am Samstag, 27. April 2019 05:43:47 UTC+2 schrieb Wolfgang Bangerth:
>>
>> On 4/26/19 2:28 PM, 'Maxi Miller' via deal.II User Group wrote: 
>> > I think the current version is the smallest version I can get. It will 
>> work 
>> > for a single thread, but not for two or more threads. 
>>
>> I tried to run this on my machine, but I don't have NOX as part of my 
>> Trilinos 
>> installation. That might have to wait for next week. 
>>
>> But I think there is still room for making the program shorter. The point 
>> is 
>> that the parts of the program that run before the error happens do not 
>> actually have to make sense. For example, could you replace your own 
>> boundary 
>> values class by ZeroFunction? Instead of assembling something real, could 
>> you 
>> just use the identity matrix on each cell? 
>>
>> Other things you don't need: Timers, convergence tables, etc. Which of 
>> the 
>> parameters you set in solve_for_NOX() are necessary? (Necessary to 
>> reproduce 
>> the bug; I know they're necessary for your application, but that's 
>> irrelevant 
>> here.) Could you replace building the residual by just putting random 
>> stuff in 
>> the vector? 
>>
>> Similarly, the computeJacobian and following functions really just check 
>> conditions, but when you run the program to illustrate the bug you're 
>> chasing, 
>> do these checks trigger? If not, remove them, then remove the calls to 
>> these 
>> functions (because they don't do anything any more), then remove the 
>> whole 
>> function. 
>>
>> I think there's still room to make this program smaller by at least a 
>> factor 
>> of 2 or 3 or 4 :-) 
>>
>> 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.
For more options, visit https://groups.google.com/d/optout.
//Nox include files
#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 
#include 
#include 

#include 
#include 
#include 

#include 
#include 

#include 
#include 
#include 

#include 

#include 
#include 

#include 
#include 

// NOX Objects
#include "NOX.H"
#include "NOX_Epetra.H"

// Trilinos Objects
#ifdef HAVE_MPI
#include "Epetra_MpiComm.h"
#else
#include "Epetra_SerialComm.h"
#endif
#include "Epetra_Vector.h"
#include "Epetra_Operator.h"
#include "Epetra_RowMatrix.h"
#include "NOX_Epetra_Interface_Required.H" // base class
#include "NOX_Epetra_Interface_Jacobian.H" // base class
#include "NOX_Epetra_Interface_Preconditioner.H" // base class
#include "Epetra_CrsMatrix.h"
#include "Epetra_Map.h"
#include "Epetra_LinearProblem.h"
#include "AztecOO.h"
#include "Teuchos_StandardCatchMacros.hpp"

// User's application specific files
//#include "problem_interface.h" // Interface file to NOX
//#include "finite_element_problem.h"

// Required for reading and writing parameter lists from xml format
// Configure Trilinos with --enable-teuchos-extended
#ifdef HAVE_TEUCHOS_EXTENDED
#include "Teuchos_XMLParameterListHelpers.hpp"
#include 
#endif

// This is C++ ...
#include 
#include 

using namespace dealii;

template 
class ProblemInterface : public NOX::Epetra::Interface::Required,
public NOX::Epetra::Interface::Preconditioner,
public NOX::Epetra::Interface::Jacobian
{
public:
ProblemInterface(std::function residual_function,
 const MPI_Comm _communicator,
   

Re: [deal.II] Re: Integrating a deal.II-specific function in a NOX-Interface for MPI-threads

['Maxi Miller' via deal.II User Group]

The functions computePreconditioner and computeJacobian are necessary, else 
I will get compilation problems. Thus I think the current version is the 
best MWE I can get. Running with one thread gives a "Test passed", using 
more threads will result in "Convergence failed".

Am Samstag, 27. April 2019 05:43:47 UTC+2 schrieb Wolfgang Bangerth:
>
> On 4/26/19 2:28 PM, 'Maxi Miller' via deal.II User Group wrote: 
> > I think the current version is the smallest version I can get. It will 
> work 
> > for a single thread, but not for two or more threads. 
>
> I tried to run this on my machine, but I don't have NOX as part of my 
> Trilinos 
> installation. That might have to wait for next week. 
>
> But I think there is still room for making the program shorter. The point 
> is 
> that the parts of the program that run before the error happens do not 
> actually have to make sense. For example, could you replace your own 
> boundary 
> values class by ZeroFunction? Instead of assembling something real, could 
> you 
> just use the identity matrix on each cell? 
>
> Other things you don't need: Timers, convergence tables, etc. Which of the 
> parameters you set in solve_for_NOX() are necessary? (Necessary to 
> reproduce 
> the bug; I know they're necessary for your application, but that's 
> irrelevant 
> here.) Could you replace building the residual by just putting random 
> stuff in 
> the vector? 
>
> Similarly, the computeJacobian and following functions really just check 
> conditions, but when you run the program to illustrate the bug you're 
> chasing, 
> do these checks trigger? If not, remove them, then remove the calls to 
> these 
> functions (because they don't do anything any more), then remove the whole 
> function. 
>
> I think there's still room to make this program smaller by at least a 
> factor 
> of 2 or 3 or 4 :-) 
>
> 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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Step 12, rewritten using ScratchData and CopyData

['Maxi Miller' via deal.II User Group]

After some tests it looks like as if the problem lies in the treatment of 
the subfaces. When refining the grid globally, everything works as 
expected, when refining adaptively, the ScratchData-Version fails. Thus I 
assume the problem lies there, but I can not see it in the code?

-- 
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] Step 12, rewritten using ScratchData and CopyData

['Maxi Miller' via deal.II User Group]

I tried to rewrite step-12 using the recently introduced classes 
ScratchData and CopyData. In the first iteration, the results are similar, 
in the second iteration the new version gives wrong results, though. I 
compared the generated matrices and right hand side vectors, and they were 
identical, but after the first mesh refinement the rewritten version 
results in wrong solutions. Thus I was wondering if I forgot to initialize 
something, or if I forgot to loop over/add some cells after the refinement?
Thanks!

-- 
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.
/* -
 *
 * Copyright (C) 2009 - 2018 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, redistribute
 * it, and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * The full text of the license can be found in the file LICENSE.md at
 * the top level directory of deal.II.
 *
 * -

 *
 * Author: Guido Kanschat, Texas A University, 2009
 */


// The first few files have already been covered in previous examples and will
// thus not be further commented on:
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
// Here the discontinuous finite elements are defined. They are used in the same
// way as all other finite elements, though -- as you have seen in previous
// tutorial programs -- there isn't much user interaction with finite element
// classes at all: they are passed to DoFHandler and
// FEValues objects, and that is about it.
#include 
// We are going to use the simplest possible solver, called Richardson
// iteration, that represents a simple defect correction. This, in combination
// with a block SSOR preconditioner (defined in precondition_block.h), that
// uses the special block matrix structure of system matrices arising from DG
// discretizations.
#include 
#include 
// We are going to use gradients as refinement indicator.
#include 

// Here come the new include files for using the MeshWorker framework. The first
// contains the class MeshWorker::DoFInfo, which provides local integrators with
// a mapping between local and global degrees of freedom. It stores the results
// of local integrals as well in its base class MeshWorker::LocalResults.
// In the second of these files, we find an object of type
// MeshWorker::IntegrationInfo, which is mostly a wrapper around a group of
// FEValues objects. The file meshworker/simple.h contains classes
// assembling locally integrated data into a global system containing only a
// single matrix. Finally, we will need the file that runs the loop over all
// mesh cells and faces.
#include 
#include 
#include 
#include 

// Like in all programs, we finish this section by including the needed C++
// headers and declaring we want to use objects in the dealii namespace without
// prefix.
#include 
#include 


namespace Step12
{
  using namespace dealii;

  // @sect3{Equation data}
  //
  // First, we define a class describing the inhomogeneous boundary data. Since
  // only its values are used, we implement value_list(), but leave all other
  // functions of Function undefined.
  template 
  class BoundaryValues : public Function
  {
  public:
BoundaryValues() = default;
virtual void value_list(const std::vector> ,
std::vector &  values,
const unsigned int component = 0) const override;
  };

  // Given the flow direction, the inflow boundary of the unit square $[0,1]^2$
  // are the right and the lower boundaries. We prescribe discontinuous boundary
  // values 1 and 0 on the x-axis and value 0 on the right boundary. The values
  // of this function on the outflow boundaries will not be used within the DG
  // scheme.
  template 
  void BoundaryValues::value_list(const std::vector> ,
   std::vector &  values,
   const unsigned int component) const
  {
(void)component;
AssertIndexRange(component, 1);
Assert(values.size() == points.size(),
   ExcDimensionMismatch(values.size(), points.size()));

for (unsigned int i = 0; i < values.size(); ++i)
  {
if 

Re: [deal.II] Re: Integrating a deal.II-specific function in a NOX-Interface for MPI-threads

['Maxi Miller' via deal.II User Group]

I think the current version is the smallest version I can get. It will work 
for a single thread, but not for two or more threads.
One observation I made was that if I not only copy the locally owned 
values, but also the locally relevant values from the deal.II-vector to the 
Epetra_Vector, the solver still does not converge, but the result is 
correct. When changing line 159 from locally_owned_elements to 
locally_relevant_elements, I get an access error, though, so there must be 
a different way to achieve that.

Am Freitag, 26. April 2019 17:00:58 UTC+2 schrieb Wolfgang Bangerth:
>
> On 4/26/19 1:59 AM, 'Maxi Miller' via deal.II User Group wrote: 
> > I tried to reduce it as much as possible while still being able to run, 
> > the result is attached. It still has ~600 lines, but I am not sure if I 
> > should remove things like the output-function? Furthermore the NOX-logic 
> > takes quite a lot of lines, and I would like to keep the program as 
> > usable as possible while removing unneccessary lines. Is that example 
> > usable? 
>
> Remove everything you can. If you get an error at one point, no output 
> will be created -- so remove the output function. If the error happens 
> before the linear system is solved, remove the code that computes the 
> entries of the matrix. Then remove the matrix object itself. If you are 
> currently solving a linear system before the error happens, try what 
> happens if you just comment out the solution step -- and if the error 
> still happens (because it probably doesn't depend on the actual values 
> of the vector), then remove the solution step and the assembly of the 
> matrix and the matrix. If you need a DoFHandler to set up the vector, 
> output the index sets and sizes you use for the processors involved and 
> build these index sets by hand instead -- then remove the DoFHandler and 
> the finite element and everything else related to this. Continue this 
> style of work until you really do have a small testcase :-) 
>
> Cheers 
>   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.
For more options, visit https://groups.google.com/d/optout.
//Nox include files
#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 
#include 
#include 

#include 
#include 
#include 

#include 
#include 

#include 
#include 
#include 

#include 

#include 
#include 

#include 
#include 

// NOX Objects
#include "NOX.H"
#include "NOX_Epetra.H"

// Trilinos Objects
#ifdef HAVE_MPI
#include "Epetra_MpiComm.h"
#else
#include "Epetra_SerialComm.h"
#endif
#include "Epetra_Vector.h"
#include "Epetra_Operator.h"
#include "Epetra_RowMatrix.h"
#include "NOX_Epetra_Interface_Required.H" // base class
#include "NOX_Epetra_Interface_Jacobian.H" // base class
#include "NOX_Epetra_Interface_Preconditioner.H" // base class
#include "Epetra_CrsMatrix.h"
#include "Epetra_Map.h"
#include "Epetra_LinearProblem.h"
#include "AztecOO.h"
#include "Teuchos_StandardCatchMacros.hpp"

// User's application specific files
//#include "problem_interface.h" // Interface file to NOX
//#include "finite_element_problem.h"

// Required for reading and writing parameter lists from xml format
// Configure Trilinos with --enable-teuchos-extended
#ifdef HAVE_TEUCHOS_EXTENDED
#include "Teuchos_XMLParameterListHelpers.hpp"
#include 
#endif

// This is C++ ...
#include 
#include 

using namespace dealii;

template 
class BoundaryValues : public dealii::Function
{
public:
BoundaryValues(const size_t n_components)
: dealii::Function(1), n_components(n_components)
{}

virtual double value (const dealii::Point   ,
  const unsigned int  component = 0) const;

virtual void vector_value(const dealii::Point , dealii::Vector ) const;
private:
const size_t n_components;
};


template 
double BoundaryValues::value (const dealii::Point ,
   const unsigned int) const
{
const double x = p[

Re: [deal.II] Re: Integrating a deal.II-specific function in a NOX-Interface for MPI-threads

['Maxi Miller' via deal.II User Group]

I tried to reduce it as much as possible while still being able to run, the 
result is attached. It still has ~600 lines, but I am not sure if I should 
remove things like the output-function? Furthermore the NOX-logic takes 
quite a lot of lines, and I would like to keep the program as usable as 
possible while removing unneccessary lines. Is that example usable?

Am Freitag, 26. April 2019 05:12:39 UTC+2 schrieb Wolfgang Bangerth:
>
> On 4/25/19 10:13 AM, 'Maxi Miller' via deal.II User Group wrote: 
> > After running some tests, I ended up reducing the problem to the 
> transfer to 
> > and from the Epetra-vectors. I have to write an interface to the model 
> > (according to the instructions), and as shown in the code above. There I 
> have 
> > Epetra-Vectors created by my interface, with a size of 
> > locally_relevant_dofs.size(), and TrilinosWrappers::MPI::Vectors. 
> Transfer 
> > from the Epetra-Vectors to the MPI::Vectors works fine, but the way back 
> does 
> > not work (The MPI::Vectors are larger than the Epetra_Vectors). Are 
> there any 
> > hints in how I still could fit the data from the MPI::Vector into the 
> > Epetra_Vector? Or should I rather ask this on the Trilinos mailing list? 
>
> Good question. I think it would probably be very useful to have small 
> testcase 
> others could look at. The program you have attached is still 1,300 lines, 
> which is far too much for anyone to understand. But since you have an idea 
> of 
> what the problem is, do you think you could come up with a small testcase 
> that 
> illustrates exactly the issue you have? It doesn't have to do anything 
> useful 
> at all -- in your case, I think all that's necessary is to create a 
> vector, 
> convert it as you describe above, and then output the sizes to show that 
> they 
> are not the same. Run this on two processors, and you should have 
> something 
> that could be done in 50 or 100 lines, and I think that might be small 
> enough 
> for someone who doesn't know your code to understand what is necessary. 
>
> I've built these sorts of testcases either from scratch in the past, or by 
> just keep removing things from a program that (i) are run after the time 
> the 
> problem happens, and (ii) then removing everything that is not necessary. 
>
> 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.
For more options, visit https://groups.google.com/d/optout.
//Nox include files
#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 
#include 
#include 

#include 
#include 
#include 

#include 
#include 

#include 
#include 
#include 

#include 

#include 
#include 

#include 
#include 

// NOX Objects
#include "NOX.H"
#include "NOX_Epetra.H"

// Trilinos Objects
#ifdef HAVE_MPI
#include "Epetra_MpiComm.h"
#else
#include "Epetra_SerialComm.h"
#endif
#include "Epetra_Vector.h"
#include "Epetra_Operator.h"
#include "Epetra_RowMatrix.h"
#include "NOX_Epetra_Interface_Required.H" // base class
#include "NOX_Epetra_Interface_Jacobian.H" // base class
#include "NOX_Epetra_Interface_Preconditioner.H" // base class
#include "Epetra_CrsMatrix.h"
#include "Epetra_Map.h"
#include "Epetra_LinearProblem.h"
#include "AztecOO.h"
#include "Teuchos_StandardCatchMacros.hpp"

// User's application specific files
//#include "problem_interface.h" // Interface file to NOX
//#include "finite_element_problem.h"

// Required for reading and writing parameter lists from xml format
// Configure Trilinos with --enable-teuchos-extended
#ifdef HAVE_TEUCHOS_EXTENDED
#include "Teuchos_XMLParameterListHelpers.hpp"
#include 
#endif

// This is C++ ...
#include 
#include 

using namespace dealii;

template 
class BoundaryValues : public dealii::Function
{
public:
BoundaryValues(const size_t n_components)
: dealii::Function(1), n_components(n_components)
{}

virtual double value (const dealii::

[deal.II] Re: Integrating a deal.II-specific function in a NOX-Interface for MPI-threads

['Maxi Miller' via deal.II User Group]

After running some tests, I ended up reducing the problem to the transfer 
to and from the Epetra-vectors. I have to write an interface to the model 
(according to the instructions), and as shown in the code above. There I 
have Epetra-Vectors created by my interface, with a size of 
locally_relevant_dofs.size(), and TrilinosWrappers::MPI::Vectors. Transfer 
from the Epetra-Vectors to the MPI::Vectors works fine, but the way back 
does not work (The MPI::Vectors are larger than the Epetra_Vectors). Are 
there any hints in how I still could fit the data from the MPI::Vector into 
the Epetra_Vector? Or should I rather ask this on the Trilinos mailing list?
Thanks!

-- 
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.
/* -
 *
 * Copyright (C) 1999 - 2018 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, redistribute
 * it, and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * The full text of the license can be found in the file LICENSE.md at
 * the top level directory of deal.II.
 *
 * -

 *
 * Author: Wolfgang Bangerth, University of Heidelberg, 1999
 */


// @sect3{Include files}

//Nox include files
#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 
#include 
#include 

#include 
#include 
#include 

#include 
#include 

#include 
#include 
#include 

#include 

#include 
#include 

#include 
#include 

// NOX Objects
#include "NOX.H"
#include "NOX_Epetra.H"

// Trilinos Objects
#ifdef HAVE_MPI
#include "Epetra_MpiComm.h"
#else
#include "Epetra_SerialComm.h"
#endif
#include "Epetra_Vector.h"
#include "Epetra_Operator.h"
#include "Epetra_RowMatrix.h"
#include "NOX_Epetra_Interface_Required.H" // base class
#include "NOX_Epetra_Interface_Jacobian.H" // base class
#include "NOX_Epetra_Interface_Preconditioner.H" // base class
#include "Epetra_CrsMatrix.h"
#include "Epetra_Map.h"
#include "Epetra_LinearProblem.h"
#include "AztecOO.h"
#include "Teuchos_StandardCatchMacros.hpp"

// User's application specific files
//#include "problem_interface.h" // Interface file to NOX
//#include "finite_element_problem.h"

// Required for reading and writing parameter lists from xml format
// Configure Trilinos with --enable-teuchos-extended
#ifdef HAVE_TEUCHOS_EXTENDED
#include "Teuchos_XMLParameterListHelpers.hpp"
#include 
#endif

// This is C++ ...
#include 
#include 

using namespace dealii;

template 
class BoundaryValues : public dealii::Function
{
public:
BoundaryValues(const size_t n_components)
: dealii::Function(1), n_components(n_components)
{}

virtual double value (const dealii::Point   ,
  const unsigned int  component = 0) const;

virtual void vector_value(const dealii::Point , dealii::Vector ) const;
private:
const size_t n_components;
};


template 
double BoundaryValues::value (const dealii::Point ,
  const unsigned int) const
{
const double x = p[0];
const double y = p[1];
return sin(M_PI * x) * sin(M_PI * y);
}

template 
void BoundaryValues::vector_value(const dealii::Point , dealii::Vector ) const
{
for(size_t i = 0; i < value.size(); ++i)
value[i] = BoundaryValues::value(p, i);
}

template 
class Solution : public Function
{
public:
Solution() : Function(1)
{

}

virtual double value(const Point , const unsigned int component) const override;
virtual Tensor<1, dim> gradient(const Point , const unsigned int component) const override;

private:
};

template 
double Solution::value(const Point , const unsigned int) const
{
const double x = p[0];
const double y = p[1];
return sin(M_PI * x) * sin(M_PI * y);
}

template
Tensor<1, dim> Solution::gradient(const Point , const unsigned int) const
{
Tensor<1, dim> return_value;
AssertThrow(dim == 2, ExcNotImplemented());

const double x = p[0];
const double y = p[1];
return_value[0] = M_PI * cos(M_PI * x) * sin(M_PI * y);
return_value[1] = M_PI * cos(M_PI * y) * sin(M_PI * x);
return return_value;
}

template 
class ProblemInterface : public 

[deal.II] Re: Error at Trilinos installation

['Maxi Miller' via deal.II User Group]

It looks as if you have to recompile Trilinos using the -fPIC-switch to 
create position independent code.

Am Donnerstag, 25. April 2019 12:10:52 UTC+2 schrieb gabrie...@koeln.de:
>
> Hey everyone,
> I have an error at compiling deal_ii_with_Trilinos=ON, which I dont 
> understand.
>
> First I installed Trilinos,then I called 
>
> cmake -DCMAKE_INSTALL_PREFIIX=/home/peters/Documents/installs 
> -DDEAL_II_WITH_TRILINOS=ON 
> -DTRILINOS_DIR=/home/peters/Documents/build_trilinos
>
> This all worked finely, but when I called make install I get the following 
> error:
>
> [ 55%] Linking CXX shared library ../lib/libdeal_II.so
> /usr/bin/ld: 
> /home/peters/Documents/build_trilinos/lib/libmuelu-adapters.a(MueLu_CreateEpetraPreconditioner.cpp.o):
>  
> relocation R_X86_64_PC32 against symbol 
> `_ZTVN7Teuchos17SerialTriDiMatrixIidEE' can not be used when making a 
> shared object; recompile with -fPIC
> /usr/bin/ld: final link failed: Bad value
> collect2: error: ld returned 1 exit status
> source/CMakeFiles/deal_II.dir/build.make:964: recipe for target 
> 'lib/libdeal_II.so.9.1.0-pre' failed
> make[2]: *** [lib/libdeal_II.so.9.1.0-pre] Error 1
> CMakeFiles/Makefile2:941: recipe for target 
> 'source/CMakeFiles/deal_II.dir/all' failed
> make[1]: *** [source/CMakeFiles/deal_II.dir/all] Error 2
> Makefile:129: recipe for target 'all' failed
> make: *** [all] Error 2
> peters@stud-04:~/programs/dealii$ m
>
>
> Can somebody tell me, what went wrong? And where do I have to add the 
> order "recompile with -fPIC".
>
> Thanks a lot and 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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Poisson test code fails for a fe-degree >= 2

['Maxi Miller' via deal.II User Group]

How exactly can I test the examples after modifying them, other than by 
running them?

Am Sonntag, 7. April 2019 23:56:49 UTC+2 schrieb Wolfgang Bangerth:
>
> On 4/7/19 1:54 PM, 'Maxi Miller' via deal.II User Group wrote: 
> > Could do that, but not until Easter (Deadlines before that), if that is 
> ok? 
>
> Any time is good for contributions! 
> Cheers 
>   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.
For more options, visit https://groups.google.com/d/optout.