Hey
I added some lines of code to the "do_integrate_difference"
function (from include/numerics/vectors.templates.h) to have
TrilinosWrappers::MPI::Vector (input vector uh) ->
TrilinosWrappers::MPI::Vector (output difference vector)
support.
See attachment.
The changes in the code are slight but not trivial,
so I'm not sure what is the best way to merge it with the current
version (explicit template instance/preprocessing statements/
or do it in another way).
If someone has an efficient idea please reply.
Best
UK
--
Dipl.-Ing. Uwe Köcher
Hamburg, Helmut Schmidt University
University of the German Federal Armed Forces
Email: [email protected]
Web: http://www.hsu-hh.de/mb-mathe
/**
* @file Integrate_Difference.tpl.cc
* @author Uwe Koecher
* @author deal.II authors
* @date 2012-03-01
*
* @brief Re-Implementation of do_integrate_difference.
*/
/* Copyright (C) 2007-2012 by the deal.II authors */
/* */
/* This file is subject to QPL and may not be distributed */
/* without copyright and license information. Please refer */
/* to the file deal.II/doc/license.html for the text and */
/* further information on this license. */
#ifndef __Integrate_Difference_tpl_cc
#define __Integrate_Difference_tpl_cc
// PROJECT includes
// MPI includes
// #include <mpi.h>
// DEAL.II includes
#include <base/utilities.h>
#include <base/function.h>
#include <base/quadrature.h>
#include <lac/vector.h>
#include <fe/mapping.h>
#include <fe/fe_values.h>
#include <lac/trilinos_vector.h>
#include <lac/constraint_matrix.h>
#include <deal.II/hp/dof_handler.h>
#include <deal.II/hp/fe_values.h>
#include <deal.II/hp/mapping_collection.h>
#include <deal.II/hp/q_collection.h>
#include <numerics/vectors.h>
using namespace dealii;
// C++ includes
#include <algorithm>
#include <cmath>
#include <limits>
#include <map>
#include <numeric>
#include <vector>
using namespace std;
#endif
namespace DTMpp {
////////////////////////////////////////////////////////////////////////////////
namespace internal {
using namespace VectorTools;
template<int dim, class InVector, class OutVector, class DH, int spacedim>
void do_integrate_difference(
const dealii::hp::MappingCollection<dim,spacedim> &mapping,
const DH &dof,
const InVector &fe_function,
const Function<spacedim> &exact_solution,
OutVector &difference,
const dealii::hp::QCollection<dim> &hp_quad,
const NormType &norm,
const Function<spacedim> *weight,
const double _exponent) {
const unsigned int n_components(dof.get_fe().n_components());
const bool fe_is_system(n_components != 1);
if (weight != NULL) {
Assert(
(weight->n_components==1) || (weight->n_components==n_components),
ExcDimensionMismatch(weight->n_components, n_components)
);
}
// difference.reinit(dof.get_tria().n_active_cells());
// Create difference vector with Epetra_Map
// (compute difference + integrate on each cell and
// save this into a difference vector);
// Each process saves the locally owned cell result.
#ifdef HAVE_MPI
Epetra_MpiComm epetra_comm(MPI_COMM_WORLD);
#else
Epetra_SerialComm epetra_comm;
#endif
const unsigned int n_global_active_cells(
dynamic_cast<const dealii::parallel::distributed::Triangulation<dim> *>(
&(dof.get_tria()))->n_global_active_cells());
const unsigned int n_locally_owned_active_cells(
dynamic_cast<const dealii::parallel::distributed::Triangulation<dim> *>(
&(dof.get_tria()))->n_locally_owned_active_cells());
Epetra_Map diff_map(n_global_active_cells, n_locally_owned_active_cells, 0,
epetra_comm);
difference.reinit(diff_map, false);
IndexSet global_index_set;
std::vector<int> global_index_vector(n_locally_owned_active_cells);
diff_map.MyGlobalElements(&(*global_index_vector.begin()));
double exponent(_exponent);
switch (norm) {
case L2_norm:
case H1_seminorm:
case H1_norm:
exponent = 2.;
break;
case L1_norm:
exponent = 1.;
break;
default:
break;
}
UpdateFlags update_flags = UpdateFlags(
update_quadrature_points | update_JxW_values);
switch (norm) {
case H1_seminorm:
case W1p_seminorm:
case W1infty_seminorm:
update_flags |= UpdateFlags(update_gradients);
if (spacedim == dim+1)
update_flags |= UpdateFlags(update_normal_vectors);
break;
case H1_norm:
case W1p_norm:
case W1infty_norm:
update_flags |= UpdateFlags(update_gradients);
if (spacedim == dim+1)
update_flags |= UpdateFlags(update_normal_vectors);
// no break!
default:
update_flags |= UpdateFlags(update_values);
break;
}
dealii::hp::FECollection<dim,spacedim> fe_collection (dof.get_fe());
dealii::hp::FEValues<dim,spacedim> x_fe_values(
mapping, fe_collection, hp_quad, update_flags);
const unsigned int max_n_q_points = hp_quad.max_n_quadrature_points();
std::vector< dealii::Vector<double> >
function_values(
max_n_q_points,
dealii::Vector<double>(n_components)
);
std::vector<std::vector<Tensor<1,spacedim> > >
function_grads(
max_n_q_points,
std::vector<Tensor<1,spacedim> >(n_components)
);
// tmp vector when we use the Function<spacedim> functions for
// scalar functions
std::vector<double> tmp_values(
max_n_q_points
);
std::vector<Tensor<1,spacedim> > tmp_gradients(
max_n_q_points
);
std::vector<double>
weight_values(
max_n_q_points
);
std::vector<dealii::Vector<double> >
weight_vectors(
max_n_q_points,
dealii::Vector<double>(n_components)
);
std::vector<dealii::Vector<double> >
psi_values(
max_n_q_points,
dealii::Vector<double>(n_components)
);
std::vector<std::vector<Tensor<1,spacedim> > >
psi_grads(
max_n_q_points,
std::vector<Tensor<1,spacedim> >(n_components)
);
std::vector<double>
psi_scalar(
max_n_q_points
);
// loop over all cells
typename DH::active_cell_iterator cell(dof.begin_active());
typename DH::active_cell_iterator endc(dof.end());
for (unsigned int index(0); cell != endc; cell++)
if (cell->is_locally_owned()) {
double diff(0);
x_fe_values.reinit(cell);
const dealii::FEValues<dim, spacedim> &fe_values(
x_fe_values.get_present_fe_values());
const unsigned int n_q_points(fe_values.n_quadrature_points);
// resize all out scratch arrays to the number of
// quadrature points we use for the present cell
function_values.resize(
n_q_points,
dealii::Vector<double>(n_components));
function_grads.resize(
n_q_points,
std::vector<Tensor<1,spacedim> >(n_components));
tmp_values.resize(
n_q_points);
tmp_gradients.resize(
n_q_points);
weight_values.resize(
n_q_points);
weight_vectors.resize(
n_q_points,
dealii::Vector<double>(n_components));
psi_values.resize(
n_q_points,
dealii::Vector<double>(n_components));
psi_grads.resize(
n_q_points,
std::vector<Tensor<1,spacedim> >(n_components));
psi_scalar.resize(
n_q_points);
if (weight != NULL) {
if (weight->n_components>1) {
weight->vector_value_list(
fe_values.get_quadrature_points(), weight_vectors);
}
else {
weight->value_list(
fe_values.get_quadrature_points(), weight_values);
for (unsigned int q(0); q < n_q_points; q++) {
weight_vectors[q] = weight_values[q];
}
}
}
else {
for (unsigned int q(0); q < n_q_points; q++) {
weight_vectors[q] = 1.;
}
}
if (update_flags & update_values) {
// first compute the exact solution
// (vectors) at the quadrature points
// try to do this as efficient as
// possible by avoiding a second
// virtual function call in case
// the function really has only
// one component
if (fe_is_system) {
exact_solution.vector_value_list(
fe_values.get_quadrature_points(), psi_values);
}
else {
exact_solution.value_list(
fe_values.get_quadrature_points(), tmp_values);
for (unsigned int q(0); q < n_q_points; q++) {
psi_values[q](0) = tmp_values[q];
}
}
// then subtract finite element
// fe_function
fe_values.get_function_values(fe_function, function_values);
for (unsigned int q(0); q < n_q_points; q++) {
psi_values[q] -= function_values[q];
}
}
// Do the same for gradients, if required
if (update_flags & update_gradients) {
// try to be a little clever
// to avoid recursive virtual
// function calls when calling
// gradient_list for functions
// that are really scalar
// functions
if (fe_is_system) {
exact_solution.vector_gradient_list(
fe_values.get_quadrature_points(), psi_grads);
}
else {
exact_solution.gradient_list(
fe_values.get_quadrature_points(), tmp_gradients);
for (unsigned int q(0); q < n_q_points; q++) {
psi_grads[q][0] = tmp_gradients[q];
}
}
// then subtract finite element function_grads. We
// need to be careful in the codimension
// one case, since there we only have tangential gradients
// in the finite element function, not the full
// gradient. This is taken care of, by subtracting the
// normal component of the gradient from the exact function.
fe_values.get_function_grads(fe_function, function_grads);
if (update_flags & update_normal_vectors) {
for (unsigned int c(0); c < n_components; c++)
for (unsigned int q(0); q < n_q_points; q++) {
psi_grads[q][c] -= (function_grads[q][c] +
(psi_grads[q][c] * // (f.n) n
fe_values.normal_vector(q)) *
fe_values.normal_vector(q));
}
}
else {
for (unsigned int c(0); c < n_components; c++)
for (unsigned int q(0); q < n_q_points; q++) {
psi_grads[q][c] -= function_grads[q][c];
}
}
}
switch (norm) {
case mean:
std::fill_n(psi_scalar.begin(), n_q_points, 0.0);
// Compute values in quadrature points
for (unsigned int c(0); c < n_components; c++)
for (unsigned int q(0); q < n_q_points; q++) {
psi_scalar[q] += psi_values[q](c) * weight_vectors[q](c);
}
// Integrate
diff = std::inner_product(psi_scalar.begin(), psi_scalar.end(),
fe_values.get_JxW_values().begin(), 0.0);
break;
case Lp_norm:
case L1_norm:
case W1p_norm:
std::fill_n(psi_scalar.begin(), n_q_points, 0.0);
// Compute values in quadrature points
for (unsigned int c(0); c < n_components; c++)
for (unsigned int q(0); q < n_q_points; q++) {
psi_scalar[q] +=
std::pow(psi_values[q](c)*psi_values[q](c), exponent/2.)
* weight_vectors[q](c);
}
// Integrate
diff = std::inner_product(psi_scalar.begin(), psi_scalar.end(),
fe_values.get_JxW_values().begin(), 0.0);
// Compute the root only,
// if no derivative values are added later
if (!(update_flags & update_gradients)) {
diff = std::pow(diff, 1./exponent);
}
break;
case L2_norm:
case H1_norm:
std::fill_n(psi_scalar.begin(), n_q_points, 0.0);
// Compute values in quadrature points
for (unsigned int c(0); c < n_components; c++)
for (unsigned int q(0); q < n_q_points; q++) {
psi_scalar[q] += psi_values[q](c)*psi_values[q](c)
* weight_vectors[q](c);
}
// Integrate
diff = std::inner_product(psi_scalar.begin(), psi_scalar.end(),
fe_values.get_JxW_values().begin(), 0.0);
// Compute the root only,
// if no derivative values are added later
if (norm == L2_norm) {
diff = std::sqrt(diff);
}
break;
case Linfty_norm:
case W1infty_norm:
std::fill_n(psi_scalar.begin(), n_q_points, 0.0);
// Compute values in quadrature points
for (unsigned int c(0); c < n_components; c++)
for (unsigned int q(0); q < n_q_points; q++) {
double newval =
std::fabs(psi_values[q](c)) * weight_vectors[q](c);
if (psi_scalar[q] < newval) {
psi_scalar[q] = newval;
}
}
// Maximum on one cell
diff = *std::max_element(psi_scalar.begin(), psi_scalar.end());
break;
case H1_seminorm:
case W1p_seminorm:
case W1infty_seminorm:
break;
default:
Assert(false, ExcNotImplemented());
break;
}
switch (norm) {
case W1p_seminorm:
case W1p_norm:
std::fill_n(psi_scalar.begin(), n_q_points, 0.0);
// Compute values in quadrature points
for (unsigned int c(0); c < n_components; c++)
for (unsigned int q(0); q < n_q_points; q++) {
psi_scalar[q] +=
std::pow(psi_grads[q][c] * psi_grads[q][c], exponent/2.)
* weight_vectors[q](c);
}
diff += std::inner_product(psi_scalar.begin(), psi_scalar.end(),
fe_values.get_JxW_values().begin(), 0.0);
diff = std::pow(diff, 1./exponent);
break;
case H1_seminorm:
case H1_norm:
// take square of integrand
std::fill_n(psi_scalar.begin(), n_q_points, 0.0);
// Compute values in quadrature points
for (unsigned int c(0); c < n_components; c++)
for (unsigned int q(0); q < n_q_points; q++) {
psi_scalar[q] +=
(psi_grads[q][c] * psi_grads[q][c]) * weight_vectors[q](c);
}
// add seminorm to L_2 norm or to zero
diff += std::inner_product(psi_scalar.begin(), psi_scalar.end(),
fe_values.get_JxW_values().begin(), 0.0);
diff = std::sqrt(diff);
break;
case W1infty_seminorm:
case W1infty_norm:
Assert(false, ExcNotImplemented());
// std::fill_n (psi_scalar.begin(), n_q_points, 0.0);
//
// // Compute values in quadrature points
// for (unsigned int c(0); c < n_components; c++)
// for (unsigned int q(0); q < n_q_points; q++) {
// double t = 0.;
// for (unsigned int d(0); d < dim; d++) {
// t = std::max(t,std::fabs(psi_grads[q][c][d])
// * weight_vectors[q](c));
// }
//
// psi_scalar[q] = std::max(psi_scalar[q],t);
// }
//
// for (unsigned int i(0); i < psi_scalar.size(); i++) {
// diff = std::max(diff, psi_scalar[i]);
// }
// break;
default:
break;
}
// append result of this cell to the end of the vector
Assert(numbers::is_finite(diff), ExcNumberNotFinite());
difference(global_index_vector[index++]) = diff;
} // if cell locally owned
// else {
// // the cell is a ghost cell or is artificial. write
// // a zero into the corresponding value of the returned vector
// difference(index) = 0;
// }
// everything is done, communicate results
difference.compress();
} // end integrate difference
} // namespace internal
} // namespace DTMpp
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
