refer to the graphical output there. The only difference between the two is
the solver and the implementation of the matrix-vector products.
-The output produced by this program is the following:
+When we run this program in 2D for quadratic elements, we get the following
+output:
@code
Cycle 0
Number of degrees of freedom: 337
@endcode
As in step-16, we see that the number of CG iterations stays constant with
-increasing number of degrees of freedom.
\ No newline at end of file
+increasing number of degrees of freedom. Not much changes if we run the
+program in three dimensions:
+@code
+Cycle 0
+Number of degrees of freedom: 517
+System matrix memory consumption: 0.09649 MBytes.
+Multigrid objects memory consumption: 0.1389 MBytes.
+Convergence in 7 CG iterations.
+
+Cycle 1
+Number of degrees of freedom: 3817
+System matrix memory consumption: 0.6333 MBytes.
+Multigrid objects memory consumption: 0.8396 MBytes.
+Convergence in 8 CG iterations.
+
+Cycle 2
+Number of degrees of freedom: 29521
+System matrix memory consumption: 4.882 MBytes.
+Multigrid objects memory consumption: 6.273 MBytes.
+Convergence in 9 CG iterations.
+
+Cycle 3
+Number of degrees of freedom: 232609
+System matrix memory consumption: 38.68 MBytes.
+Multigrid objects memory consumption: 49.26 MBytes.
+Convergence in 10 CG iterations.
+
+Cycle 4
+Number of degrees of freedom: 1847617
+System matrix memory consumption: 308.4 MBytes.
+Multigrid objects memory consumption: 391.5 MBytes.
+Convergence in 11 CG iterations.
+@endcode
+
+There is a slight increase in the number of CG iterations, but also for the
+largest size with almost two million unknowns, the iteration count is still
+very modest.
\ No newline at end of file
LaplaceOperator<dim,number>&
operator = (const Tensor<2,dim> &tensor);
- number transformation[dim][dim];
+ number transformation[dim*(dim+1)/2];
};
template<int dim,typename number>
temp_result[d] = result[d];
for (unsigned int d=0; d<dim; ++d)
{
- number output = 0;
- for (unsigned int e=0; e<dim; ++e)
- output += transformation[d][e] * temp_result[e];
+ number output = transformation[d]*temp_result[d];
+ if (dim == 2)
+ output += transformation[2]*temp_result[1-d];
+ else if (dim == 3)
+ {
+ if (d==0)
+ output += transformation[3]*temp_result[1] + transformation[4]*temp_result[2];
+ else if (d==1)
+ output += transformation[3]*temp_result[0] + transformation[5]*temp_result[2];
+ else
+ output += transformation[4]*temp_result[0] + transformation[5]*temp_result[1];
+ }
result[d] = output;
}
}
template <int dim, typename number>
LaplaceOperator<dim,number>&
-LaplaceOperator<dim,number>::operator=(const Tensor<2,dim> &tensor)
+LaplaceOperator<dim,number>::operator=(const Tensor<2,dim> &tensor)
{
- for (unsigned int d=0;d<dim;++d)
- for (unsigned int e=0;e<dim;++e)
- transformation[d][e] = tensor[d][e];
+ if (dim == 2)
+ {
+ transformation[0] = tensor[0][0];
+ transformation[1] = tensor[1][1];
+ transformation[2] = tensor[0][1];
+ Assert (std::fabs(tensor[1][0]-tensor[0][1])<1e-15,
+ ExcInternalError());
+ }
+ else if (dim == 3)
+ {
+ transformation[0] = tensor[0][0];
+ transformation[1] = tensor[1][1];
+ transformation[2] = tensor[2][2];
+ transformation[3] = tensor[0][1];
+ transformation[4] = tensor[0][2];
+ transformation[5] = tensor[1][2];
+ Assert (std::fabs(tensor[1][0]-tensor[0][1])<1e-15,
+ ExcInternalError());
+ Assert (std::fabs(tensor[2][0]-tensor[0][2])<1e-15,
+ ExcInternalError());
+ Assert (std::fabs(tensor[2][1]-tensor[1][2])<1e-15,
+ ExcInternalError());
+ }
return *this;
}
template <int dim>
void LaplaceProblem<dim>::run ()
{
- for (unsigned int cycle=0; cycle<6; ++cycle)
+ for (unsigned int cycle=0; cycle<8-dim; ++cycle)
{
std::cout << "Cycle " << cycle << std::endl;
GridGenerator::hyper_ball(triangulation);
static const HyperBallBoundary<dim> boundary;
triangulation.set_boundary (0, boundary);
- triangulation.refine_global (0);
+ triangulation.refine_global (3-dim);
}
triangulation.refine_global (1);
setup_system ();