From 846311e3b77a953d15f3035649fc735c8709d899 Mon Sep 17 00:00:00 2001 From: guido Date: Tue, 30 Aug 2005 09:48:05 +0000 Subject: [PATCH] move instatiations to vectors.templates.h; thanks to Dima Sorkin git-svn-id: https://svn.dealii.org/trunk@11339 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/numerics/vectors.h | 16 + .../include/numerics/vectors.templates.h | 1648 ++++++++++++++ deal.II/deal.II/source/numerics/vectors.cc | 1964 +---------------- .../source/numerics/vectors.instance.h | 103 + 4 files changed, 1782 insertions(+), 1949 deletions(-) create mode 100644 deal.II/deal.II/include/numerics/vectors.templates.h create mode 100644 deal.II/deal.II/source/numerics/vectors.instance.h diff --git a/deal.II/deal.II/include/numerics/vectors.h b/deal.II/deal.II/include/numerics/vectors.h index 22447c2012..2b37acf61a 100644 --- a/deal.II/deal.II/include/numerics/vectors.h +++ b/deal.II/deal.II/include/numerics/vectors.h @@ -292,6 +292,7 @@ class ConstraintMatrix; * if access to an object describing the exact form of the boundary is needed, the * pointer stored within the triangulation object is accessed. * + * @ref Instantiations: some (Vector, Vector, BlockVector, BlockVector, see also individual functions) * @author Wolfgang Bangerth, Ralf Hartmann, Guido Kanschat, 1998, 1999, 2000, 2001 */ class VectorTools @@ -432,6 +433,15 @@ class VectorTools * space @p dof_2 afterwards, to * make the result continuous * again. + * + * @ref Instantiations: some (the + * standard vector types + * mentioned in the class + * documentation are + * instantiated, but only the + * same vector for InVector and + * OutVector. Other combinations + * must be instantiated by hand.) */ template static void interpolate (const DoFHandler &dof_1, @@ -830,6 +840,12 @@ class VectorTools * * See the general documentation of this * class for more information. + * + * @ref Instantiations: some + * (InVectors as in the + * documentation of the class, + * OutVector only Vector + * and Vector) */ template static void integrate_difference (const Mapping &mapping, diff --git a/deal.II/deal.II/include/numerics/vectors.templates.h b/deal.II/deal.II/include/numerics/vectors.templates.h new file mode 100644 index 0000000000..4f07b685c2 --- /dev/null +++ b/deal.II/deal.II/include/numerics/vectors.templates.h @@ -0,0 +1,1648 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2005 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 _deal2__vectors_templates_h +#define _deal2__vectors_templates_h + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +template +inline double sqr_point (const Tensor<1,dim> &p) +{ + return p * p; +} + + +template +void VectorTools::interpolate (const Mapping &mapping, + const DoFHandler &dof, + const Function &function, + VECTOR &vec) +{ + Assert (dof.get_fe().n_components() == function.n_components, + ExcComponentMismatch()); + + const FiniteElement &fe = dof.get_fe(); + const unsigned int n_components = fe.n_components(); + const bool fe_is_system = (n_components != 1); + + typename DoFHandler::active_cell_iterator cell = dof.begin_active(), + endc = dof.end(); + + // For FESystems many of the + // unit_support_points will + // appear multiply, as a point + // may be unit_support_point + // for several of the components + // of the system. + // The following is rather + // complicated as it is + // avoided to evaluate + // the vectorfunction multiply at + // the same point on a cell. + const std::vector > & + unit_support_points = fe.get_unit_support_points(); + Assert (unit_support_points.size() != 0, + ExcNonInterpolatingFE()); + + // Find the support points + // on a cell that + // are multiply mentioned in + // @p{unit_support_points}. + // Mark the first representative + // of each multiply mentioned + // support point by appending its + // dof index to @p{dofs_of_rep_points}. + // Each multiple point gets to know + // the dof index of its representative + // point by the @p{dof_to_rep_dof_table}. + + // the following vector collects all dofs i, + // 0<=i dofs_of_rep_points; + // the following table converts a dof i + // to the rep index. + std::vector dof_to_rep_index_table; + unsigned int n_rep_points=0; + for (unsigned int i=0; i0; --j) + if (unit_support_points[i] + == unit_support_points[dofs_of_rep_points[j-1]]) + { + dof_to_rep_index_table.push_back(j-1); + representative=false; + break; + } + + if (representative) + { + // rep_index=dofs_of_rep_points.size() + dof_to_rep_index_table.push_back(dofs_of_rep_points.size()); + // dofs_of_rep_points[rep_index]=i + dofs_of_rep_points.push_back(i); + ++n_rep_points; + } + } + Assert(dofs_of_rep_points.size()==n_rep_points, ExcInternalError()); + Assert(dof_to_rep_index_table.size()==fe.dofs_per_cell, ExcInternalError()); + + std::vector dofs_on_cell (fe.dofs_per_cell); + std::vector > rep_points (n_rep_points); + + // get space for the values of the + // function at the rep support points. + // + // have two versions, one for system fe + // and one for scalar ones, to take the + // more efficient one respectively + std::vector function_values_scalar (n_rep_points); + std::vector > function_values_system (n_rep_points, + Vector(fe.n_components())); + + // Make a quadrature rule from support points + // to feed it into FEValues + Quadrature support_quadrature(unit_support_points); + + // Transformed support points are computed by + // FEValues + FEValues fe_values (mapping, fe, support_quadrature, update_q_points); + + for (; cell!=endc; ++cell) + { + // for each cell: + // get location of finite element + // support_points + fe_values.reinit(cell); + const std::vector >& support_points = + fe_values.get_quadrature_points(); + + // pick out the representative + // support points + for (unsigned int j=0; jget_dof_indices (dofs_on_cell); + + + if (fe_is_system) + { + // get function values at + // these points. Here: get + // all components + function.vector_value_list (rep_points, function_values_system); + // distribute the function + // values to the global + // vector + for (unsigned int i=0; i +void VectorTools::interpolate (const DoFHandler &dof, + const Function &function, + VECTOR &vec) +{ + Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); + static const MappingQ1 mapping; + interpolate(mapping, dof, function, vec); +} + + + + +template +void +VectorTools::interpolate (const DoFHandler &dof_1, + const DoFHandler &dof_2, + const FullMatrix &transfer, + const InVector &data_1, + OutVector &data_2) +{ + Vector cell_data_1(dof_1.get_fe().dofs_per_cell); + Vector cell_data_2(dof_2.get_fe().dofs_per_cell); + + std::vector touch_count (dof_2.n_dofs(), 0); + std::vector local_dof_indices (dof_2.get_fe().dofs_per_cell); + + typename DoFHandler::active_cell_iterator h = dof_1.begin_active(); + typename DoFHandler::active_cell_iterator l = dof_2.begin_active(); + const typename DoFHandler::cell_iterator endh = dof_1.end(); + + for(; h != endh; ++h, ++l) + { + h->get_dof_values(data_1, cell_data_1); + transfer.vmult(cell_data_2, cell_data_1); + + l->get_dof_indices (local_dof_indices); + + // distribute cell vector + for (unsigned int j=0; j &, + const DoFHandler<1> &, + const ConstraintMatrix &, + const Quadrature<1> &, + const Function<1> &, + Vector &, + const bool , + const Quadrature<0> &, + const bool ) +{ + // this function should easily be implemented + // using the template below. However some + // changes have to be made since faces don't + // exist in 1D. Maybe integrate the creation of + // zero boundary values into the + // project_boundary_values function? + Assert (false, ExcNotImplemented()); +} + + +#endif + + +template +void VectorTools::project (const Mapping &mapping, + const DoFHandler &dof, + const ConstraintMatrix &constraints, + const Quadrature &quadrature, + const Function &function, + Vector &vec, + const bool enforce_zero_boundary, + const Quadrature &q_boundary, + const bool project_to_boundary_first) +{ + Assert (dof.get_fe().n_components() == function.n_components, + ExcInvalidFE()); + + const FiniteElement &fe = dof.get_fe(); + + // make up boundary values + std::map boundary_values; + + if (enforce_zero_boundary == true) + // no need to project boundary + // values, but enforce + // homogeneous boundary values + // anyway + { + // loop over all boundary faces + // to get all dof indices of + // dofs on the boundary. note + // that in 3d there are cases + // where a face is not at the + // boundary, yet one of its + // lines is, and we should + // consider the degrees of + // freedom on it as boundary + // nodes. likewise, in 2d and + // 3d there are cases where a + // cell is only at the boundary + // by one vertex. nevertheless, + // since we do not support + // boundaries with dimension + // less or equal to dim-2, each + // such boundary dof is also + // found from some other face + // that is actually wholly on + // the boundary, not only by + // one line or one vertex + typename DoFHandler::active_face_iterator face = dof.begin_active_face(), + endf = dof.end_face(); + std::vector face_dof_indices (fe.dofs_per_face); + for (; face!=endf; ++face) + if (face->at_boundary()) + { + face->get_dof_indices (face_dof_indices); + for (unsigned int i=0; i::type boundary_functions; + for (unsigned char c=0; c<255; ++c) + boundary_functions[c] = &function; + project_boundary_values (dof, boundary_functions, q_boundary, + boundary_values); + }; + + + // set up mass matrix and right hand side + vec.reinit (dof.n_dofs()); + SparsityPattern sparsity(dof.n_dofs(), + dof.n_dofs(), + dof.max_couplings_between_dofs()); + DoFTools::make_sparsity_pattern (dof, sparsity); + constraints.condense (sparsity); + + SparseMatrix mass_matrix (sparsity); + Vector tmp (mass_matrix.n()); + + MatrixCreator::create_mass_matrix (mapping, dof, quadrature, mass_matrix); + + VectorTools::create_right_hand_side (mapping, dof, quadrature, function, tmp); + + constraints.condense (mass_matrix); + constraints.condense (tmp); + if (boundary_values.size() != 0) + MatrixTools::apply_boundary_values (boundary_values, + mass_matrix, vec, tmp, + true); + + SolverControl control(1000,1e-16); + PrimitiveVectorMemory<> memory; + SolverCG<> cg(control,memory); + + PreconditionSSOR<> prec; + prec.initialize(mass_matrix, 1.2); + // solve + cg.solve (mass_matrix, vec, tmp, prec); + + // distribute solution + constraints.distribute (vec); +} + + +template +void VectorTools::project (const DoFHandler &dof, + const ConstraintMatrix &constraints, + const Quadrature &quadrature, + const Function &function, + Vector &vec, + const bool enforce_zero_boundary, + const Quadrature &q_boundary, + const bool project_to_boundary_first) +{ + Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); + static const MappingQ1 mapping; + project(mapping, dof, constraints, quadrature, function, vec, + enforce_zero_boundary, q_boundary, project_to_boundary_first); +} + + + + +template +void VectorTools::create_right_hand_side (const Mapping &mapping, + const DoFHandler &dof_handler, + const Quadrature &quadrature, + const Function &rhs_function, + Vector &rhs_vector) +{ + const FiniteElement &fe = dof_handler.get_fe(); + Assert (fe.n_components() == rhs_function.n_components, + ExcComponentMismatch()); + Assert (rhs_vector.size() == dof_handler.n_dofs(), + ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs())); + rhs_vector = 0; + + UpdateFlags update_flags = UpdateFlags(update_values | + update_q_points | + update_JxW_values); + FEValues fe_values (mapping, fe, quadrature, update_flags); + + const unsigned int dofs_per_cell = fe_values.dofs_per_cell, + n_q_points = fe_values.n_quadrature_points, + n_components = fe.n_components(); + + std::vector dofs (dofs_per_cell); + Vector cell_vector (dofs_per_cell); + + typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), + endc = dof_handler.end(); + + if (n_components==1) + { + std::vector rhs_values(n_q_points); + + for (; cell!=endc; ++cell) + { + fe_values.reinit(cell); + + const std::vector &weights = fe_values.get_JxW_values (); + rhs_function.value_list (fe_values.get_quadrature_points(), rhs_values); + + cell_vector = 0; + for (unsigned int point=0; pointget_dof_indices (dofs); + + for (unsigned int i=0; i > rhs_values(n_q_points, Vector(n_components)); + + // Use the faster code if the FiniteElement is primitive + if (fe.is_primitive ()) + { + for (; cell!=endc; ++cell) + { + fe_values.reinit(cell); + + const std::vector &weights = fe_values.get_JxW_values (); + rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values); + + cell_vector = 0; + for (unsigned int point=0; pointget_dof_indices (dofs); + + for (unsigned int i=0; i &weights = fe_values.get_JxW_values (); + rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values); + + cell_vector = 0; + for (unsigned int point=0; pointget_dof_indices (dofs); + + for (unsigned int i=0; i +void VectorTools::create_right_hand_side (const DoFHandler &dof_handler, + const Quadrature &quadrature, + const Function &rhs_function, + Vector &rhs_vector) +{ + Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); + static const MappingQ1 mapping; + create_right_hand_side(mapping, dof_handler, quadrature, + rhs_function, rhs_vector); +} + + + +#if deal_II_dimension != 1 + +template +void +VectorTools::create_boundary_right_hand_side (const Mapping &mapping, + const DoFHandler &dof_handler, + const Quadrature &quadrature, + const Function &rhs_function, + Vector &rhs_vector, + const std::set &boundary_indicators) +{ + const FiniteElement &fe = dof_handler.get_fe(); + Assert (fe.n_components() == rhs_function.n_components, + ExcComponentMismatch()); + Assert (rhs_vector.size() == dof_handler.n_dofs(), + ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs())); + + rhs_vector = 0; + + UpdateFlags update_flags = UpdateFlags(update_values | + update_q_points | + update_JxW_values); + FEFaceValues fe_values (mapping, fe, quadrature, update_flags); + + const unsigned int dofs_per_cell = fe_values.dofs_per_cell, + n_q_points = fe_values.n_quadrature_points, + n_components = fe.n_components(); + + std::vector dofs (dofs_per_cell); + Vector cell_vector (dofs_per_cell); + + typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), + endc = dof_handler.end(); + + if (n_components==1) + { + std::vector rhs_values(n_q_points); + + for (; cell!=endc; ++cell) + for (unsigned int face=0; face::faces_per_cell; ++face) + if (cell->face(face)->at_boundary () && + (boundary_indicators.find (cell->face(face)->boundary_indicator()) + != + boundary_indicators.end())) + { + fe_values.reinit(cell, face); + + const std::vector &weights = fe_values.get_JxW_values (); + rhs_function.value_list (fe_values.get_quadrature_points(), rhs_values); + + cell_vector = 0; + for (unsigned int point=0; pointget_dof_indices (dofs); + + for (unsigned int i=0; i > rhs_values(n_q_points, Vector(n_components)); + + for (; cell!=endc; ++cell) + for (unsigned int face=0; face::faces_per_cell; ++face) + if (cell->face(face)->at_boundary () && + (boundary_indicators.find (cell->face(face)->boundary_indicator()) + != + boundary_indicators.end())) + { + fe_values.reinit(cell, face); + + const std::vector &weights = fe_values.get_JxW_values (); + rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values); + + cell_vector = 0; + + // Use the faster code if the FiniteElement is primitive + if (fe.is_primitive ()) + { + for (unsigned int point=0; pointget_dof_indices (dofs); + + for (unsigned int i=0; i &, + const DoFHandler<1> &, + const Quadrature<0> &, + const Function<1> &, + Vector &, + const std::set &) +{ + Assert (false, ExcImpossibleInDim(1)); +} + +#endif + +template +void +VectorTools::create_boundary_right_hand_side (const DoFHandler &dof_handler, + const Quadrature &quadrature, + const Function &rhs_function, + Vector &rhs_vector, + const std::set &boundary_indicators) +{ + Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); + static const MappingQ1 mapping; + create_boundary_right_hand_side(mapping, dof_handler, quadrature, + rhs_function, rhs_vector, + boundary_indicators); +} + + + +#if deal_II_dimension == 1 + +void +VectorTools::interpolate_boundary_values (const Mapping<1> &, + const DoFHandler<1> &dof, + const unsigned char boundary_component, + const Function<1> &boundary_function, + std::map &boundary_values, + const std::vector &component_mask_) +{ + Assert (boundary_component != 255, + ExcInvalidBoundaryIndicator()); + + const FiniteElement<1> &fe = dof.get_fe(); + Assert (fe.n_components() == boundary_function.n_components, + ExcComponentMismatch()); + + // set the component mask to either + // the original value or a vector + // of @p{true}s + const std::vector component_mask ((component_mask_.size() == 0) ? + std::vector (fe.n_components(), true) : + component_mask_); + Assert (std::count(component_mask.begin(), component_mask.end(), true) > 0, + ExcComponentMismatch()); + + // check whether boundary values at + // the left or right boundary of + // the line are + // requested. @p{direction} denotes + // the neighboring direction in + // which we seek the boundary, + // i.e. 0 is left boundary and 1 is + // right. + const unsigned int direction = boundary_component; + Assert (direction < 2, ExcInvalidBoundaryIndicator()); + + // first find the outermost active + // cell by first traversing the coarse + // grid to its end and then going + // to the children + DoFHandler<1>::cell_iterator outermost_cell = dof.begin(0); + while (outermost_cell->neighbor(direction).state() == IteratorState::valid) + outermost_cell = outermost_cell->neighbor(direction); + + while (outermost_cell->has_children()) + outermost_cell = outermost_cell->child(direction); + + // now set the value of the + // outermost degree of + // freedom. setting also + // creates the entry in the map + // if it did not exist + // beforehand + // + // save some time by requesting + // values only once for each point, + // irrespective of the number of + // components of the function + Vector function_values (fe.n_components()); + if (fe.n_components() == 1) + function_values(0) + = boundary_function.value (outermost_cell->vertex(direction)); + else + boundary_function.vector_value (outermost_cell->vertex(direction), + function_values); + + for (unsigned int i=0; ivertex_dof_index(direction,i)] + = function_values(fe.face_system_to_component_index(i).first); +} + + + +void +VectorTools::interpolate_boundary_values (const Mapping<1> &mapping, + const DoFHandler<1> &dof, + const FunctionMap<1>::type &function_map, + std::map &boundary_values, + const std::vector &component_mask) +{ + for (FunctionMap<1>::type::const_iterator i=function_map.begin(); + i!=function_map.end(); ++i) + interpolate_boundary_values (mapping, dof, i->first, *i->second, + boundary_values, component_mask); +} + + +#endif + + +template +void +VectorTools:: +interpolate_boundary_values (const Mapping &mapping, + const DoFHandler &dof, + const typename FunctionMap::type &function_map, + std::map &boundary_values, + const std::vector &component_mask_) +{ + // if for whatever reason we were + // passed an empty map, return + // immediately + if (function_map.size() == 0) + return; + + Assert (function_map.find(255) == function_map.end(), + ExcInvalidBoundaryIndicator()); + + const FiniteElement &fe = dof.get_fe(); + const unsigned int n_components = fe.n_components(); + const bool fe_is_system = (n_components != 1); + + for (typename FunctionMap::type::const_iterator i=function_map.begin(); + i!=function_map.end(); ++i) + Assert (n_components == i->second->n_components, + ExcInvalidFE()); + + // set the component mask to either + // the original value or a vector + // of @p{true}s + const std::vector component_mask ((component_mask_.size() == 0) ? + std::vector (fe.n_components(), true) : + component_mask_); + Assert (std::count(component_mask.begin(), component_mask.end(), true) > 0, + ExcComponentMismatch()); + + // field to store the indices + std::vector face_dofs (fe.dofs_per_face, + DoFHandler::invalid_dof_index); + std::vector > dof_locations (face_dofs.size(), Point()); + + // array to store the values of + // the boundary function at the + // boundary points. have to arrays + // for scalar and vector functions + // to use the more efficient one + // respectively + std::vector dof_values_scalar (fe.dofs_per_face); + std::vector > dof_values_system (fe.dofs_per_face, + Vector(fe.n_components())); + + // next generate a quadrature rule + // on the face from the unit + // support points. this wil be used + // to obtain the quadrature points + // on the real cell's face + std::vector > + unit_support_points = fe.get_unit_face_support_points(); + + // check whether there are support + // points on the face. if not, then + // we should try a more clever + // way. the idea is that a finite + // element may not offer support + // points for all its shape + // functions, but maybe only + // some. if it offers support + // points for the components we are + // interested in in this function, + // then that's fine. if not, the + // function we call in the finite + // element will raise an + // exception. the support points + // for the other shape functions + // are left uninitialized (well, + // initialized by the default + // constructor), since we don't + // need them anyway. + if (unit_support_points.size() == 0) + { + unit_support_points.resize (fe.dofs_per_face); + for (unsigned int i=0; i aux_quad (unit_support_points); + FEFaceValues fe_values (mapping, fe, aux_quad, update_q_points); + + typename DoFHandler::active_cell_iterator cell = dof.begin_active(), + endc = dof.end(); + for (; cell!=endc; ++cell) + for (unsigned int face_no = 0; face_no < GeometryInfo::faces_per_cell; + ++face_no) + { + typename DoFHandler::face_iterator face = cell->face(face_no); + const unsigned char boundary_component = face->boundary_indicator(); + if (function_map.find(boundary_component) != function_map.end()) + // face is of the right component + { + // get indices, physical location and + // boundary values of dofs on this + // face + face->get_dof_indices (face_dofs); + fe_values.reinit(cell, face_no); + const std::vector > &dof_locations = fe_values.get_quadrature_points (); + + if (fe_is_system) + { + function_map.find(boundary_component)->second + ->vector_value_list (dof_locations, dof_values_system); + + // enter those dofs + // into the list that + // match the + // component + // signature. avoid + // the usual + // complication that + // we can't just use + // *_system_to_component_index + // for non-primitive + // FEs + for (unsigned int i=0; isecond->value_list (dof_locations, + dof_values_scalar, + 0); + + // enter into list + + for (unsigned int i=0; i +void +VectorTools::interpolate_boundary_values (const Mapping &mapping, + const DoFHandler &dof, + const unsigned char boundary_component, + const Function &boundary_function, + std::map &boundary_values, + const std::vector &component_mask) +{ + typename FunctionMap::type function_map; + function_map[boundary_component] = &boundary_function; + interpolate_boundary_values (mapping, dof, function_map, boundary_values, + component_mask); +} + + + +template +void +VectorTools::interpolate_boundary_values (const DoFHandler &dof, + const unsigned char boundary_component, + const Function &boundary_function, + std::map &boundary_values, + const std::vector &component_mask) +{ + Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); + static const MappingQ1 mapping; + interpolate_boundary_values(mapping, dof, boundary_component, + boundary_function, boundary_values, component_mask); +} + + + +template +void +VectorTools::interpolate_boundary_values (const DoFHandler &dof, + const typename FunctionMap::type &function_map, + std::map &boundary_values, + const std::vector &component_mask) +{ + Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); + static const MappingQ1 mapping; + interpolate_boundary_values(mapping, dof, function_map, + boundary_values, component_mask); +} + + +#if deal_II_dimension == 1 + +void +VectorTools::project_boundary_values (const Mapping<1> &mapping, + const DoFHandler<1> &dof, + const FunctionMap<1>::type &boundary_functions, + const Quadrature<0> &, + std::map &boundary_values) +{ + // projection in 1d is equivalent + // to interpolation + interpolate_boundary_values (mapping, dof, boundary_functions, + boundary_values, std::vector()); +} + +#endif + + +template +void +VectorTools::project_boundary_values (const Mapping &mapping, + const DoFHandler &dof, + const typename FunctionMap::type &boundary_functions, + const Quadrature &q, + std::map &boundary_values) +{ +//TODO:[?] In VectorTools::project_boundary_values, no condensation of sparsity +// structures, matrices and right hand sides or distribution of +// solution vectors is performed. This is ok for dim<3 because then +// there are no constrained nodes on the boundary, but is not +// acceptable for higher dimensions. Fix this. + + Assert (dof.get_fe().n_components() == boundary_functions.begin()->second->n_components, + ExcComponentMismatch()); + + std::vector dof_to_boundary_mapping; + std::set selected_boundary_components; + for (typename FunctionMap::type::const_iterator i=boundary_functions.begin(); + i!=boundary_functions.end(); ++i) + selected_boundary_components.insert (i->first); + + DoFTools::map_dof_to_boundary_indices (dof, selected_boundary_components, + dof_to_boundary_mapping); + + // set up sparsity structure + SparsityPattern sparsity(dof.n_boundary_dofs(boundary_functions), + dof.max_couplings_between_boundary_dofs()); + DoFTools::make_boundary_sparsity_pattern (dof, + boundary_functions, + dof_to_boundary_mapping, + sparsity); + + // note: for three or more dimensions, there + // may be constrained nodes on the boundary + // in this case the boundary mass matrix has + // to be condensed and the solution is to + // be distributed afterwards, which is not + // yet implemented. The reason for this is + // that we cannot simply use the @p{condense} + // family of functions, since the matrices + // and vectors do not use the global + // numbering but rather the boundary + // numbering, i.e. the condense function + // needs to use another indirection. There + // should be not many technical problems, + // but it needs to be implemented + if (dim<3) + sparsity.compress(); + else + Assert (false, ExcNotImplemented()); + + + // make mass matrix and right hand side + SparseMatrix mass_matrix(sparsity); + Vector rhs(sparsity.n_rows()); + + + MatrixCreator::create_boundary_mass_matrix (mapping, dof, q, + mass_matrix, boundary_functions, + rhs, dof_to_boundary_mapping); + + // same thing as above: if dim>=3 we need + // to consider constraints + Assert (dim<3, ExcNotImplemented()); + + + Vector boundary_projection (rhs.size()); + + SolverControl control(1000, 1e-16); + PrimitiveVectorMemory<> memory; + SolverCG<> cg(control,memory); + + PreconditionSSOR<> prec; + prec.initialize(mass_matrix, 1.2); + // solve + cg.solve (mass_matrix, boundary_projection, rhs, prec); + + // fill in boundary values + for (unsigned int i=0; i::invalid_dof_index) + // this dof is on one of the + // interesting boundary parts + // + // remember: @p{i} is the global dof + // number, @p{dof_to_boundary_mapping[i]} + // is the number on the boundary and + // thus in the solution vector + boundary_values[i] = boundary_projection(dof_to_boundary_mapping[i]); +} + + +template +void +VectorTools::project_boundary_values (const DoFHandler &dof, + const typename FunctionMap::type &boundary_functions, + const Quadrature &q, + std::map &boundary_values) +{ + Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); + static const MappingQ1 mapping; + project_boundary_values(mapping, dof, boundary_functions, q, boundary_values); +} + + + +template +void +VectorTools::integrate_difference (const Mapping &mapping, + const DoFHandler &dof, + const InVector &fe_function, + const Function &exact_solution, + OutVector &difference, + const Quadrature &q, + const NormType &norm, + const Function *weight, + const double exponent_1) +{ + // we mark the "exponent" parameter + // to this function "const" since + // it is strictly incoming, but we + // need to set it to something + // different later on, if + // necessary, so have a read-write + // version of it: + double exponent = exponent_1; + + const unsigned int n_q_points = q.n_quadrature_points; + const FiniteElement &fe = dof.get_fe(); + const unsigned int n_components = fe.n_components(); + const bool fe_is_system = (n_components != 1); + + if (weight!=0) + { + Assert ((weight->n_components==1) || (weight->n_components==n_components), + ExcDimensionMismatch(weight->n_components, n_components)); + } + + difference.reinit (dof.get_tria().n_active_cells()); + + 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_q_points | + update_JxW_values); + switch (norm) + { + case H1_seminorm: + case W1p_seminorm: + case W1infty_seminorm: + update_flags |= UpdateFlags (update_gradients); + break; + case H1_norm: + case W1p_norm: + case W1infty_norm: + update_flags |= UpdateFlags (update_gradients); + // no break! + default: + update_flags |= UpdateFlags (update_values); + break; + } + + FEValues fe_values(mapping, fe, q, update_flags); + + std::vector< Vector > function_values (n_q_points, + Vector(n_components)); + std::vector > > function_grads (n_q_points, + std::vector >(n_components)); + std::vector weight_values (n_q_points); + std::vector > weight_vectors (n_q_points, + Vector(n_components)); + + std::vector > psi_values (n_q_points, + Vector(n_components)); + std::vector > > psi_grads (n_q_points, + std::vector >(n_components)); + std::vector psi_scalar (n_q_points); + // tmp vector when we use the + // Function functions for + // scalar functions + std::vector tmp_values (fe_values.n_quadrature_points); + std::vector > tmp_gradients (fe_values.n_quadrature_points); + + // loop over all cells + typename DoFHandler::active_cell_iterator cell = dof.begin_active(), + endc = dof.end(); + for (unsigned int index=0; cell != endc; ++cell, ++index) + { + double diff=0; + // initialize for this cell + fe_values.reinit (cell); + + if (weight!=0) + { + 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 k=0;k +void +VectorTools::integrate_difference (const DoFHandler &dof, + const InVector &fe_function, + const Function &exact_solution, + OutVector &difference, + const Quadrature &q, + const NormType &norm, + const Function *weight, + const double exponent) +{ + Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); + static const MappingQ1 mapping; + integrate_difference(mapping, dof, fe_function, exact_solution, + difference, q, norm, weight, exponent); +} + + + +template +void +VectorTools::point_difference (const DoFHandler &dof, + const InVector &fe_function, + const Function &exact_function, + Vector &difference, + const Point &point) +{ + static const MappingQ1 mapping; + const FiniteElement& fe = dof.get_fe(); + + Assert(difference.size() == fe.n_components(), + ExcDimensionMismatch(difference.size(), fe.n_components())); + + // first find the cell in which this point + // is, initialize a quadrature rule with + // it, and then a FEValues object + const typename DoFHandler::active_cell_iterator + cell = GridTools::find_active_cell_around_point (dof, point); + + const Point unit_point + = mapping.transform_real_to_unit_cell(cell, point); + Assert (GeometryInfo::is_inside_unit_cell (unit_point), + ExcInternalError()); + + const Quadrature quadrature (std::vector > (1, unit_point), + std::vector (1, 1.)); + FEValues fe_values(mapping, fe, quadrature, update_values); + fe_values.reinit(cell); + + // then use this to get at the values of + // the given fe_function at this point + std::vector > u_value(1, Vector (fe.n_components())); + fe_values.get_function_values(fe_function, u_value); + + if (fe.n_components() == 1) + difference(0) = exact_function.value(point); + else + exact_function.vector_value(point, difference); + + for (unsigned int i=0; i +double +VectorTools::compute_mean_value (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &quadrature, + const InVector &v, + const unsigned int component) +{ + Assert (component < dof.get_fe().n_components(), + ExcIndexRange(component, 0, dof.get_fe().n_components())); + + FEValues fe(mapping, dof.get_fe(), quadrature, + UpdateFlags(update_JxW_values + | update_values)); + + typename DoFHandler::active_cell_iterator c; + std::vector > values(quadrature.n_quadrature_points, + Vector (dof.get_fe().n_components())); + + double mean = 0.; + double area = 0.; + // Compute mean value + for (c = dof.begin_active(); c != dof.end(); ++c) + { + fe.reinit (c); + fe.get_function_values(v, values); + for (unsigned int k=0; k< quadrature.n_quadrature_points; ++k) + { + mean += fe.JxW(k) * values[k](component); + area += fe.JxW(k); + }; + }; + + return (mean/area); +} + + +template +double +VectorTools::compute_mean_value (const DoFHandler &dof, + const Quadrature &quadrature, + const InVector &v, + const unsigned int component) +{ + Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); + static const MappingQ1 mapping; + return compute_mean_value(mapping, dof, quadrature, v, component); +} + +#endif diff --git a/deal.II/deal.II/source/numerics/vectors.cc b/deal.II/deal.II/source/numerics/vectors.cc index 2bbad617ac..3dbc4d1779 100644 --- a/deal.II/deal.II/source/numerics/vectors.cc +++ b/deal.II/deal.II/source/numerics/vectors.cc @@ -11,1705 +11,25 @@ // //--------------------------------------------------------------------------- +#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include - -//TODO:[GK] Move templates containing vector arguments to vectors.templates.h - - - -template -inline double sqr_point (const Tensor<1,dim> &p) -{ - return p * p; -} - - - - -template -void VectorTools::interpolate (const Mapping &mapping, - const DoFHandler &dof, - const Function &function, - VECTOR &vec) -{ - Assert (dof.get_fe().n_components() == function.n_components, - ExcComponentMismatch()); - - const FiniteElement &fe = dof.get_fe(); - const unsigned int n_components = fe.n_components(); - const bool fe_is_system = (n_components != 1); - - typename DoFHandler::active_cell_iterator cell = dof.begin_active(), - endc = dof.end(); - - // For FESystems many of the - // unit_support_points will - // appear multiply, as a point - // may be unit_support_point - // for several of the components - // of the system. - // The following is rather - // complicated as it is - // avoided to evaluate - // the vectorfunction multiply at - // the same point on a cell. - const std::vector > & - unit_support_points = fe.get_unit_support_points(); - Assert (unit_support_points.size() != 0, - ExcNonInterpolatingFE()); - - // Find the support points - // on a cell that - // are multiply mentioned in - // @p{unit_support_points}. - // Mark the first representative - // of each multiply mentioned - // support point by appending its - // dof index to @p{dofs_of_rep_points}. - // Each multiple point gets to know - // the dof index of its representative - // point by the @p{dof_to_rep_dof_table}. - - // the following vector collects all dofs i, - // 0<=i dofs_of_rep_points; - // the following table converts a dof i - // to the rep index. - std::vector dof_to_rep_index_table; - unsigned int n_rep_points=0; - for (unsigned int i=0; i0; --j) - if (unit_support_points[i] - == unit_support_points[dofs_of_rep_points[j-1]]) - { - dof_to_rep_index_table.push_back(j-1); - representative=false; - break; - } - - if (representative) - { - // rep_index=dofs_of_rep_points.size() - dof_to_rep_index_table.push_back(dofs_of_rep_points.size()); - // dofs_of_rep_points[rep_index]=i - dofs_of_rep_points.push_back(i); - ++n_rep_points; - } - } - Assert(dofs_of_rep_points.size()==n_rep_points, ExcInternalError()); - Assert(dof_to_rep_index_table.size()==fe.dofs_per_cell, ExcInternalError()); - - std::vector dofs_on_cell (fe.dofs_per_cell); - std::vector > rep_points (n_rep_points); - - // get space for the values of the - // function at the rep support points. - // - // have two versions, one for system fe - // and one for scalar ones, to take the - // more efficient one respectively - std::vector function_values_scalar (n_rep_points); - std::vector > function_values_system (n_rep_points, - Vector(fe.n_components())); - - // Make a quadrature rule from support points - // to feed it into FEValues - Quadrature support_quadrature(unit_support_points); - - // Transformed support points are computed by - // FEValues - FEValues fe_values (mapping, fe, support_quadrature, update_q_points); - - for (; cell!=endc; ++cell) - { - // for each cell: - // get location of finite element - // support_points - fe_values.reinit(cell); - const std::vector >& support_points = - fe_values.get_quadrature_points(); - - // pick out the representative - // support points - for (unsigned int j=0; jget_dof_indices (dofs_on_cell); - - - if (fe_is_system) - { - // get function values at - // these points. Here: get - // all components - function.vector_value_list (rep_points, function_values_system); - // distribute the function - // values to the global - // vector - for (unsigned int i=0; i -void VectorTools::interpolate (const DoFHandler &dof, - const Function &function, - VECTOR &vec) -{ - Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); - static const MappingQ1 mapping; - interpolate(mapping, dof, function, vec); -} - - - - -template -void -VectorTools::interpolate (const DoFHandler &dof_1, - const DoFHandler &dof_2, - const FullMatrix &transfer, - const InVector &data_1, - OutVector &data_2) -{ - Vector cell_data_1(dof_1.get_fe().dofs_per_cell); - Vector cell_data_2(dof_2.get_fe().dofs_per_cell); - - std::vector touch_count (dof_2.n_dofs(), 0); - std::vector local_dof_indices (dof_2.get_fe().dofs_per_cell); - - typename DoFHandler::active_cell_iterator h = dof_1.begin_active(); - typename DoFHandler::active_cell_iterator l = dof_2.begin_active(); - const typename DoFHandler::cell_iterator endh = dof_1.end(); - - for(; h != endh; ++h, ++l) - { - h->get_dof_values(data_1, cell_data_1); - transfer.vmult(cell_data_2, cell_data_1); - - l->get_dof_indices (local_dof_indices); - - // distribute cell vector - for (unsigned int j=0; j &, - const DoFHandler<1> &, - const ConstraintMatrix &, - const Quadrature<1> &, - const Function<1> &, - Vector &, - const bool , - const Quadrature<0> &, - const bool ) -{ - // this function should easily be implemented - // using the template below. However some - // changes have to be made since faces don't - // exist in 1D. Maybe integrate the creation of - // zero boundary values into the - // project_boundary_values function? - Assert (false, ExcNotImplemented()); -} - - -#endif - - -template -void VectorTools::project (const Mapping &mapping, - const DoFHandler &dof, - const ConstraintMatrix &constraints, - const Quadrature &quadrature, - const Function &function, - Vector &vec, - const bool enforce_zero_boundary, - const Quadrature &q_boundary, - const bool project_to_boundary_first) -{ - Assert (dof.get_fe().n_components() == function.n_components, - ExcInvalidFE()); - - const FiniteElement &fe = dof.get_fe(); - - // make up boundary values - std::map boundary_values; - - if (enforce_zero_boundary == true) - // no need to project boundary - // values, but enforce - // homogeneous boundary values - // anyway - { - // loop over all boundary faces - // to get all dof indices of - // dofs on the boundary. note - // that in 3d there are cases - // where a face is not at the - // boundary, yet one of its - // lines is, and we should - // consider the degrees of - // freedom on it as boundary - // nodes. likewise, in 2d and - // 3d there are cases where a - // cell is only at the boundary - // by one vertex. nevertheless, - // since we do not support - // boundaries with dimension - // less or equal to dim-2, each - // such boundary dof is also - // found from some other face - // that is actually wholly on - // the boundary, not only by - // one line or one vertex - typename DoFHandler::active_face_iterator face = dof.begin_active_face(), - endf = dof.end_face(); - std::vector face_dof_indices (fe.dofs_per_face); - for (; face!=endf; ++face) - if (face->at_boundary()) - { - face->get_dof_indices (face_dof_indices); - for (unsigned int i=0; i::type boundary_functions; - for (unsigned char c=0; c<255; ++c) - boundary_functions[c] = &function; - project_boundary_values (dof, boundary_functions, q_boundary, - boundary_values); - }; - - - // set up mass matrix and right hand side - vec.reinit (dof.n_dofs()); - SparsityPattern sparsity(dof.n_dofs(), - dof.n_dofs(), - dof.max_couplings_between_dofs()); - DoFTools::make_sparsity_pattern (dof, sparsity); - constraints.condense (sparsity); - - SparseMatrix mass_matrix (sparsity); - Vector tmp (mass_matrix.n()); - - MatrixCreator::create_mass_matrix (mapping, dof, quadrature, mass_matrix); - - VectorTools::create_right_hand_side (mapping, dof, quadrature, function, tmp); - - constraints.condense (mass_matrix); - constraints.condense (tmp); - if (boundary_values.size() != 0) - MatrixTools::apply_boundary_values (boundary_values, - mass_matrix, vec, tmp, - true); - - SolverControl control(1000,1e-16); - PrimitiveVectorMemory<> memory; - SolverCG<> cg(control,memory); - - PreconditionSSOR<> prec; - prec.initialize(mass_matrix, 1.2); - // solve - cg.solve (mass_matrix, vec, tmp, prec); - - // distribute solution - constraints.distribute (vec); -} - - -template -void VectorTools::project (const DoFHandler &dof, - const ConstraintMatrix &constraints, - const Quadrature &quadrature, - const Function &function, - Vector &vec, - const bool enforce_zero_boundary, - const Quadrature &q_boundary, - const bool project_to_boundary_first) -{ - Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); - static const MappingQ1 mapping; - project(mapping, dof, constraints, quadrature, function, vec, - enforce_zero_boundary, q_boundary, project_to_boundary_first); -} - - - - -template -void VectorTools::create_right_hand_side (const Mapping &mapping, - const DoFHandler &dof_handler, - const Quadrature &quadrature, - const Function &rhs_function, - Vector &rhs_vector) -{ - const FiniteElement &fe = dof_handler.get_fe(); - Assert (fe.n_components() == rhs_function.n_components, - ExcComponentMismatch()); - Assert (rhs_vector.size() == dof_handler.n_dofs(), - ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs())); - rhs_vector = 0; - - UpdateFlags update_flags = UpdateFlags(update_values | - update_q_points | - update_JxW_values); - FEValues fe_values (mapping, fe, quadrature, update_flags); - - const unsigned int dofs_per_cell = fe_values.dofs_per_cell, - n_q_points = fe_values.n_quadrature_points, - n_components = fe.n_components(); - - std::vector dofs (dofs_per_cell); - Vector cell_vector (dofs_per_cell); - - typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), - endc = dof_handler.end(); - - if (n_components==1) - { - std::vector rhs_values(n_q_points); - - for (; cell!=endc; ++cell) - { - fe_values.reinit(cell); - - const std::vector &weights = fe_values.get_JxW_values (); - rhs_function.value_list (fe_values.get_quadrature_points(), rhs_values); - - cell_vector = 0; - for (unsigned int point=0; pointget_dof_indices (dofs); - - for (unsigned int i=0; i > rhs_values(n_q_points, Vector(n_components)); - - // Use the faster code if the FiniteElement is primitive - if (fe.is_primitive ()) - { - for (; cell!=endc; ++cell) - { - fe_values.reinit(cell); - - const std::vector &weights = fe_values.get_JxW_values (); - rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values); - - cell_vector = 0; - for (unsigned int point=0; pointget_dof_indices (dofs); - - for (unsigned int i=0; i &weights = fe_values.get_JxW_values (); - rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values); - - cell_vector = 0; - for (unsigned int point=0; pointget_dof_indices (dofs); - - for (unsigned int i=0; i -void VectorTools::create_right_hand_side (const DoFHandler &dof_handler, - const Quadrature &quadrature, - const Function &rhs_function, - Vector &rhs_vector) -{ - Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); - static const MappingQ1 mapping; - create_right_hand_side(mapping, dof_handler, quadrature, - rhs_function, rhs_vector); -} - - - -#if deal_II_dimension != 1 - -template -void -VectorTools::create_boundary_right_hand_side (const Mapping &mapping, - const DoFHandler &dof_handler, - const Quadrature &quadrature, - const Function &rhs_function, - Vector &rhs_vector, - const std::set &boundary_indicators) -{ - const FiniteElement &fe = dof_handler.get_fe(); - Assert (fe.n_components() == rhs_function.n_components, - ExcComponentMismatch()); - Assert (rhs_vector.size() == dof_handler.n_dofs(), - ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs())); - - rhs_vector = 0; - - UpdateFlags update_flags = UpdateFlags(update_values | - update_q_points | - update_JxW_values); - FEFaceValues fe_values (mapping, fe, quadrature, update_flags); - - const unsigned int dofs_per_cell = fe_values.dofs_per_cell, - n_q_points = fe_values.n_quadrature_points, - n_components = fe.n_components(); - - std::vector dofs (dofs_per_cell); - Vector cell_vector (dofs_per_cell); - - typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), - endc = dof_handler.end(); - - if (n_components==1) - { - std::vector rhs_values(n_q_points); - - for (; cell!=endc; ++cell) - for (unsigned int face=0; face::faces_per_cell; ++face) - if (cell->face(face)->at_boundary () && - (boundary_indicators.find (cell->face(face)->boundary_indicator()) - != - boundary_indicators.end())) - { - fe_values.reinit(cell, face); - - const std::vector &weights = fe_values.get_JxW_values (); - rhs_function.value_list (fe_values.get_quadrature_points(), rhs_values); - - cell_vector = 0; - for (unsigned int point=0; pointget_dof_indices (dofs); - - for (unsigned int i=0; i > rhs_values(n_q_points, Vector(n_components)); - - for (; cell!=endc; ++cell) - for (unsigned int face=0; face::faces_per_cell; ++face) - if (cell->face(face)->at_boundary () && - (boundary_indicators.find (cell->face(face)->boundary_indicator()) - != - boundary_indicators.end())) - { - fe_values.reinit(cell, face); - - const std::vector &weights = fe_values.get_JxW_values (); - rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values); - - cell_vector = 0; - - // Use the faster code if the FiniteElement is primitive - if (fe.is_primitive ()) - { - for (unsigned int point=0; pointget_dof_indices (dofs); - - for (unsigned int i=0; i &, - const DoFHandler<1> &, - const Quadrature<0> &, - const Function<1> &, - Vector &, - const std::set &) -{ - Assert (false, ExcImpossibleInDim(1)); -} - -#endif - -template -void -VectorTools::create_boundary_right_hand_side (const DoFHandler &dof_handler, - const Quadrature &quadrature, - const Function &rhs_function, - Vector &rhs_vector, - const std::set &boundary_indicators) -{ - Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); - static const MappingQ1 mapping; - create_boundary_right_hand_side(mapping, dof_handler, quadrature, - rhs_function, rhs_vector, - boundary_indicators); -} - - - -#if deal_II_dimension == 1 - -void -VectorTools::interpolate_boundary_values (const Mapping<1> &, - const DoFHandler<1> &dof, - const unsigned char boundary_component, - const Function<1> &boundary_function, - std::map &boundary_values, - const std::vector &component_mask_) -{ - Assert (boundary_component != 255, - ExcInvalidBoundaryIndicator()); - - const FiniteElement<1> &fe = dof.get_fe(); - Assert (fe.n_components() == boundary_function.n_components, - ExcComponentMismatch()); - - // set the component mask to either - // the original value or a vector - // of @p{true}s - const std::vector component_mask ((component_mask_.size() == 0) ? - std::vector (fe.n_components(), true) : - component_mask_); - Assert (std::count(component_mask.begin(), component_mask.end(), true) > 0, - ExcComponentMismatch()); - - // check whether boundary values at - // the left or right boundary of - // the line are - // requested. @p{direction} denotes - // the neighboring direction in - // which we seek the boundary, - // i.e. 0 is left boundary and 1 is - // right. - const unsigned int direction = boundary_component; - Assert (direction < 2, ExcInvalidBoundaryIndicator()); - - // first find the outermost active - // cell by first traversing the coarse - // grid to its end and then going - // to the children - DoFHandler<1>::cell_iterator outermost_cell = dof.begin(0); - while (outermost_cell->neighbor(direction).state() == IteratorState::valid) - outermost_cell = outermost_cell->neighbor(direction); - - while (outermost_cell->has_children()) - outermost_cell = outermost_cell->child(direction); - - // now set the value of the - // outermost degree of - // freedom. setting also - // creates the entry in the map - // if it did not exist - // beforehand - // - // save some time by requesting - // values only once for each point, - // irrespective of the number of - // components of the function - Vector function_values (fe.n_components()); - if (fe.n_components() == 1) - function_values(0) - = boundary_function.value (outermost_cell->vertex(direction)); - else - boundary_function.vector_value (outermost_cell->vertex(direction), - function_values); - - for (unsigned int i=0; ivertex_dof_index(direction,i)] - = function_values(fe.face_system_to_component_index(i).first); -} - - - -void -VectorTools::interpolate_boundary_values (const Mapping<1> &mapping, - const DoFHandler<1> &dof, - const FunctionMap<1>::type &function_map, - std::map &boundary_values, - const std::vector &component_mask) -{ - for (FunctionMap<1>::type::const_iterator i=function_map.begin(); - i!=function_map.end(); ++i) - interpolate_boundary_values (mapping, dof, i->first, *i->second, - boundary_values, component_mask); -} - - -#endif - - -template -void -VectorTools:: -interpolate_boundary_values (const Mapping &mapping, - const DoFHandler &dof, - const typename FunctionMap::type &function_map, - std::map &boundary_values, - const std::vector &component_mask_) -{ - // if for whatever reason we were - // passed an empty map, return - // immediately - if (function_map.size() == 0) - return; - - Assert (function_map.find(255) == function_map.end(), - ExcInvalidBoundaryIndicator()); - - const FiniteElement &fe = dof.get_fe(); - const unsigned int n_components = fe.n_components(); - const bool fe_is_system = (n_components != 1); - - for (typename FunctionMap::type::const_iterator i=function_map.begin(); - i!=function_map.end(); ++i) - Assert (n_components == i->second->n_components, - ExcInvalidFE()); - - // set the component mask to either - // the original value or a vector - // of @p{true}s - const std::vector component_mask ((component_mask_.size() == 0) ? - std::vector (fe.n_components(), true) : - component_mask_); - Assert (std::count(component_mask.begin(), component_mask.end(), true) > 0, - ExcComponentMismatch()); - - // field to store the indices - std::vector face_dofs (fe.dofs_per_face, - DoFHandler::invalid_dof_index); - std::vector > dof_locations (face_dofs.size(), Point()); - - // array to store the values of - // the boundary function at the - // boundary points. have to arrays - // for scalar and vector functions - // to use the more efficient one - // respectively - std::vector dof_values_scalar (fe.dofs_per_face); - std::vector > dof_values_system (fe.dofs_per_face, - Vector(fe.n_components())); - - // next generate a quadrature rule - // on the face from the unit - // support points. this wil be used - // to obtain the quadrature points - // on the real cell's face - std::vector > - unit_support_points = fe.get_unit_face_support_points(); - - // check whether there are support - // points on the face. if not, then - // we should try a more clever - // way. the idea is that a finite - // element may not offer support - // points for all its shape - // functions, but maybe only - // some. if it offers support - // points for the components we are - // interested in in this function, - // then that's fine. if not, the - // function we call in the finite - // element will raise an - // exception. the support points - // for the other shape functions - // are left uninitialized (well, - // initialized by the default - // constructor), since we don't - // need them anyway. - if (unit_support_points.size() == 0) - { - unit_support_points.resize (fe.dofs_per_face); - for (unsigned int i=0; i aux_quad (unit_support_points); - FEFaceValues fe_values (mapping, fe, aux_quad, update_q_points); - - typename DoFHandler::active_cell_iterator cell = dof.begin_active(), - endc = dof.end(); - for (; cell!=endc; ++cell) - for (unsigned int face_no = 0; face_no < GeometryInfo::faces_per_cell; - ++face_no) - { - typename DoFHandler::face_iterator face = cell->face(face_no); - const unsigned char boundary_component = face->boundary_indicator(); - if (function_map.find(boundary_component) != function_map.end()) - // face is of the right component - { - // get indices, physical location and - // boundary values of dofs on this - // face - face->get_dof_indices (face_dofs); - fe_values.reinit(cell, face_no); - const std::vector > &dof_locations = fe_values.get_quadrature_points (); - - if (fe_is_system) - { - function_map.find(boundary_component)->second - ->vector_value_list (dof_locations, dof_values_system); - - // enter those dofs - // into the list that - // match the - // component - // signature. avoid - // the usual - // complication that - // we can't just use - // *_system_to_component_index - // for non-primitive - // FEs - for (unsigned int i=0; isecond->value_list (dof_locations, - dof_values_scalar, - 0); - - // enter into list - - for (unsigned int i=0; i -void -VectorTools::interpolate_boundary_values (const Mapping &mapping, - const DoFHandler &dof, - const unsigned char boundary_component, - const Function &boundary_function, - std::map &boundary_values, - const std::vector &component_mask) -{ - typename FunctionMap::type function_map; - function_map[boundary_component] = &boundary_function; - interpolate_boundary_values (mapping, dof, function_map, boundary_values, - component_mask); -} - - - -template -void -VectorTools::interpolate_boundary_values (const DoFHandler &dof, - const unsigned char boundary_component, - const Function &boundary_function, - std::map &boundary_values, - const std::vector &component_mask) -{ - Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); - static const MappingQ1 mapping; - interpolate_boundary_values(mapping, dof, boundary_component, - boundary_function, boundary_values, component_mask); -} - - - -template -void -VectorTools::interpolate_boundary_values (const DoFHandler &dof, - const typename FunctionMap::type &function_map, - std::map &boundary_values, - const std::vector &component_mask) -{ - Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); - static const MappingQ1 mapping; - interpolate_boundary_values(mapping, dof, function_map, - boundary_values, component_mask); -} - - -#if deal_II_dimension == 1 - -void -VectorTools::project_boundary_values (const Mapping<1> &mapping, - const DoFHandler<1> &dof, - const FunctionMap<1>::type &boundary_functions, - const Quadrature<0> &, - std::map &boundary_values) -{ - // projection in 1d is equivalent - // to interpolation - interpolate_boundary_values (mapping, dof, boundary_functions, - boundary_values, std::vector()); -} - -#endif - - -template -void -VectorTools::project_boundary_values (const Mapping &mapping, - const DoFHandler &dof, - const typename FunctionMap::type &boundary_functions, - const Quadrature &q, - std::map &boundary_values) -{ -//TODO:[?] In VectorTools::project_boundary_values, no condensation of sparsity -// structures, matrices and right hand sides or distribution of -// solution vectors is performed. This is ok for dim<3 because then -// there are no constrained nodes on the boundary, but is not -// acceptable for higher dimensions. Fix this. - - Assert (dof.get_fe().n_components() == boundary_functions.begin()->second->n_components, - ExcComponentMismatch()); - - std::vector dof_to_boundary_mapping; - std::set selected_boundary_components; - for (typename FunctionMap::type::const_iterator i=boundary_functions.begin(); - i!=boundary_functions.end(); ++i) - selected_boundary_components.insert (i->first); - - DoFTools::map_dof_to_boundary_indices (dof, selected_boundary_components, - dof_to_boundary_mapping); - - // set up sparsity structure - SparsityPattern sparsity(dof.n_boundary_dofs(boundary_functions), - dof.max_couplings_between_boundary_dofs()); - DoFTools::make_boundary_sparsity_pattern (dof, - boundary_functions, - dof_to_boundary_mapping, - sparsity); - - // note: for three or more dimensions, there - // may be constrained nodes on the boundary - // in this case the boundary mass matrix has - // to be condensed and the solution is to - // be distributed afterwards, which is not - // yet implemented. The reason for this is - // that we cannot simply use the @p{condense} - // family of functions, since the matrices - // and vectors do not use the global - // numbering but rather the boundary - // numbering, i.e. the condense function - // needs to use another indirection. There - // should be not many technical problems, - // but it needs to be implemented - if (dim<3) - sparsity.compress(); - else - Assert (false, ExcNotImplemented()); - - - // make mass matrix and right hand side - SparseMatrix mass_matrix(sparsity); - Vector rhs(sparsity.n_rows()); - - - MatrixCreator::create_boundary_mass_matrix (mapping, dof, q, - mass_matrix, boundary_functions, - rhs, dof_to_boundary_mapping); - - // same thing as above: if dim>=3 we need - // to consider constraints - Assert (dim<3, ExcNotImplemented()); - - - Vector boundary_projection (rhs.size()); - - SolverControl control(1000, 1e-16); - PrimitiveVectorMemory<> memory; - SolverCG<> cg(control,memory); - - PreconditionSSOR<> prec; - prec.initialize(mass_matrix, 1.2); - // solve - cg.solve (mass_matrix, boundary_projection, rhs, prec); - - // fill in boundary values - for (unsigned int i=0; i::invalid_dof_index) - // this dof is on one of the - // interesting boundary parts - // - // remember: @p{i} is the global dof - // number, @p{dof_to_boundary_mapping[i]} - // is the number on the boundary and - // thus in the solution vector - boundary_values[i] = boundary_projection(dof_to_boundary_mapping[i]); -} - - -template -void -VectorTools::project_boundary_values (const DoFHandler &dof, - const typename FunctionMap::type &boundary_functions, - const Quadrature &q, - std::map &boundary_values) -{ - Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); - static const MappingQ1 mapping; - project_boundary_values(mapping, dof, boundary_functions, q, boundary_values); -} - - - -template -void -VectorTools::integrate_difference (const Mapping &mapping, - const DoFHandler &dof, - const InVector &fe_function, - const Function &exact_solution, - OutVector &difference, - const Quadrature &q, - const NormType &norm, - const Function *weight, - const double exponent_1) -{ - // we mark the "exponent" parameter - // to this function "const" since - // it is strictly incoming, but we - // need to set it to something - // different later on, if - // necessary, so have a read-write - // version of it: - double exponent = exponent_1; - - const unsigned int n_q_points = q.n_quadrature_points; - const FiniteElement &fe = dof.get_fe(); - const unsigned int n_components = fe.n_components(); - const bool fe_is_system = (n_components != 1); - - if (weight!=0) - { - Assert ((weight->n_components==1) || (weight->n_components==n_components), - ExcDimensionMismatch(weight->n_components, n_components)); - } - - difference.reinit (dof.get_tria().n_active_cells()); - - 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_q_points | - update_JxW_values); - switch (norm) - { - case H1_seminorm: - case W1p_seminorm: - case W1infty_seminorm: - update_flags |= UpdateFlags (update_gradients); - break; - case H1_norm: - case W1p_norm: - case W1infty_norm: - update_flags |= UpdateFlags (update_gradients); - // no break! - default: - update_flags |= UpdateFlags (update_values); - break; - } - - FEValues fe_values(mapping, fe, q, update_flags); - - std::vector< Vector > function_values (n_q_points, - Vector(n_components)); - std::vector > > function_grads (n_q_points, - std::vector >(n_components)); - std::vector weight_values (n_q_points); - std::vector > weight_vectors (n_q_points, - Vector(n_components)); - - std::vector > psi_values (n_q_points, - Vector(n_components)); - std::vector > > psi_grads (n_q_points, - std::vector >(n_components)); - std::vector psi_scalar (n_q_points); - // tmp vector when we use the - // Function functions for - // scalar functions - std::vector tmp_values (fe_values.n_quadrature_points); - std::vector > tmp_gradients (fe_values.n_quadrature_points); - - // loop over all cells - typename DoFHandler::active_cell_iterator cell = dof.begin_active(), - endc = dof.end(); - for (unsigned int index=0; cell != endc; ++cell, ++index) - { - double diff=0; - // initialize for this cell - fe_values.reinit (cell); - - if (weight!=0) - { - 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 k=0;k -void -VectorTools::integrate_difference (const DoFHandler &dof, - const InVector &fe_function, - const Function &exact_solution, - OutVector &difference, - const Quadrature &q, - const NormType &norm, - const Function *weight, - const double exponent) -{ - Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); - static const MappingQ1 mapping; - integrate_difference(mapping, dof, fe_function, exact_solution, - difference, q, norm, weight, exponent); -} - - - -template -void -VectorTools::point_difference (const DoFHandler &dof, - const InVector &fe_function, - const Function &exact_function, - Vector &difference, - const Point &point) -{ - static const MappingQ1 mapping; - const FiniteElement& fe = dof.get_fe(); - - Assert(difference.size() == fe.n_components(), - ExcDimensionMismatch(difference.size(), fe.n_components())); - - // first find the cell in which this point - // is, initialize a quadrature rule with - // it, and then a FEValues object - const typename DoFHandler::active_cell_iterator - cell = GridTools::find_active_cell_around_point (dof, point); - - const Point unit_point - = mapping.transform_real_to_unit_cell(cell, point); - Assert (GeometryInfo::is_inside_unit_cell (unit_point), - ExcInternalError()); - - const Quadrature quadrature (std::vector > (1, unit_point), - std::vector (1, 1.)); - FEValues fe_values(mapping, fe, quadrature, update_values); - fe_values.reinit(cell); - - // then use this to get at the values of - // the given fe_function at this point - std::vector > u_value(1, Vector (fe.n_components())); - fe_values.get_function_values(fe_function, u_value); - - if (fe.n_components() == 1) - difference(0) = exact_function.value(point); - else - exact_function.vector_value(point, difference); - - for (unsigned int i=0; i -double -VectorTools::compute_mean_value (const Mapping &mapping, - const DoFHandler &dof, - const Quadrature &quadrature, - const InVector &v, - const unsigned int component) -{ - Assert (component < dof.get_fe().n_components(), - ExcIndexRange(component, 0, dof.get_fe().n_components())); - - FEValues fe(mapping, dof.get_fe(), quadrature, - UpdateFlags(update_JxW_values - | update_values)); +// explicit instantiations - typename DoFHandler::active_cell_iterator c; - std::vector > values(quadrature.n_quadrature_points, - Vector (dof.get_fe().n_components())); - - double mean = 0.; - double area = 0.; - // Compute mean value - for (c = dof.begin_active(); c != dof.end(); ++c) - { - fe.reinit (c); - fe.get_function_values(v, values); - for (unsigned int k=0; k< quadrature.n_quadrature_points; ++k) - { - mean += fe.JxW(k) * values[k](component); - area += fe.JxW(k); - }; - }; - - return (mean/area); -} +#define VEC Vector +#include "vectors.instance.h" +#undef VEC +#define VEC Vector +#include "vectors.instance.h" +#undef VEC -template -double -VectorTools::compute_mean_value (const DoFHandler &dof, - const Quadrature &quadrature, - const InVector &v, - const unsigned int component) -{ - Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); - static const MappingQ1 mapping; - return compute_mean_value(mapping, dof, quadrature, v, component); -} +#define VEC BlockVector +#include "vectors.instance.h" +#undef VEC - -// explicit instantiations -template -void VectorTools::interpolate -(const Mapping&, - const DoFHandler&, - const Function&, - Vector&); -template -void VectorTools::interpolate -(const DoFHandler&, - const Function&, - Vector&); -template -void VectorTools::interpolate -(const Mapping&, - const DoFHandler&, - const Function&, - Vector&); -template -void VectorTools::interpolate -(const DoFHandler&, - const Function&, - Vector&); -template -void VectorTools::interpolate -(const Mapping&, - const DoFHandler&, - const Function&, - BlockVector&); -template -void VectorTools::interpolate -(const DoFHandler&, - const Function&, - BlockVector&); -template -void VectorTools::interpolate -(const Mapping&, - const DoFHandler&, - const Function&, - BlockVector&); -template -void VectorTools::interpolate -(const DoFHandler&, - const Function&, - BlockVector&); - -template -void VectorTools::interpolate -(const DoFHandler&, - const DoFHandler&, - const FullMatrix&, - const Vector&, - Vector&); -template -void VectorTools::interpolate -(const DoFHandler&, - const DoFHandler&, - const FullMatrix&, - const BlockVector&, - BlockVector&); +#define VEC BlockVector +#include "vectors.instance.h" +#undef VEC template void VectorTools::project @@ -1721,6 +41,7 @@ void VectorTools::project const bool, const Quadrature &, const bool); + template void VectorTools::create_right_hand_side (const Mapping &, @@ -1762,204 +83,6 @@ void VectorTools::interpolate_boundary_values ( const Function &, std::map &, const std::vector &); -template -void VectorTools::integrate_difference ( - const Mapping &, - const DoFHandler &, - const Vector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); -template -void VectorTools::integrate_difference -(const DoFHandler &, - const Vector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); -template -void VectorTools::integrate_difference -(const Mapping &, - const DoFHandler &, - const Vector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); -template -void VectorTools::integrate_difference -(const DoFHandler &, - const Vector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); -template -void VectorTools::integrate_difference -(const Mapping &, - const DoFHandler &, - const Vector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); -template -void VectorTools::integrate_difference -(const DoFHandler &, - const Vector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); -template -void VectorTools::integrate_difference -(const Mapping &, - const DoFHandler &, - const Vector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); -template -void VectorTools::integrate_difference -(const DoFHandler &, - const Vector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); -template -void VectorTools::integrate_difference -(const Mapping &, - const DoFHandler &, - const BlockVector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); -template -void VectorTools::integrate_difference -(const DoFHandler &, - const BlockVector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); -template -void VectorTools::integrate_difference -(const Mapping &, - const DoFHandler &, - const BlockVector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); -template -void VectorTools::integrate_difference -(const DoFHandler &, - const BlockVector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); -template -void VectorTools::integrate_difference -(const Mapping &, - const DoFHandler &, - const BlockVector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); -template -void VectorTools::integrate_difference -(const DoFHandler &, - const BlockVector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); -template -void VectorTools::integrate_difference -(const Mapping &, - const DoFHandler &, - const BlockVector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); -template -void VectorTools::integrate_difference -(const DoFHandler &, - const BlockVector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); - -template -void VectorTools::point_difference ( - const DoFHandler&, - const Vector&, - const Function&, - Vector&, - const Point&); -template -void VectorTools::point_difference ( - const DoFHandler&, - const Vector&, - const Function&, - Vector&, - const Point&); -template -void VectorTools::point_difference ( - const DoFHandler&, - const BlockVector&, - const Function&, - Vector&, - const Point&); -template -void VectorTools::point_difference ( - const DoFHandler&, - const BlockVector&, - const Function&, - Vector&, - const Point&); - #if deal_II_dimension != 1 template @@ -1978,63 +101,6 @@ void VectorTools::project_boundary_values const Quadrature&, std::map &); -template -double VectorTools::compute_mean_value -(const Mapping&, - const DoFHandler&, - const Quadrature&, - const Vector&, - const unsigned int); -template -double VectorTools::compute_mean_value -(const DoFHandler&, - const Quadrature&, - const Vector&, - const unsigned int); - -template -double VectorTools::compute_mean_value -(const Mapping&, - const DoFHandler&, - const Quadrature&, - const BlockVector&, - const unsigned int); -template -double VectorTools::compute_mean_value -(const DoFHandler&, - const Quadrature&, - const BlockVector&, - const unsigned int); - -template -double VectorTools::compute_mean_value -(const Mapping&, - const DoFHandler&, - const Quadrature&, - const Vector&, - const unsigned int); -template -double VectorTools::compute_mean_value -(const DoFHandler&, - const Quadrature&, - const Vector&, - const unsigned int); - -template -double VectorTools::compute_mean_value -(const Mapping&, - const DoFHandler&, - const Quadrature&, - const BlockVector&, - const unsigned int); -template -double VectorTools::compute_mean_value -(const DoFHandler&, - const Quadrature&, - const BlockVector&, - const unsigned int); - - // the following two functions are not derived from a template in 1d // and thus need no explicit instantiation diff --git a/deal.II/deal.II/source/numerics/vectors.instance.h b/deal.II/deal.II/source/numerics/vectors.instance.h new file mode 100644 index 0000000000..a2d1464b22 --- /dev/null +++ b/deal.II/deal.II/source/numerics/vectors.instance.h @@ -0,0 +1,103 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2005 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. +// +//--------------------------------------------------------------------------- + +// Instantiations of functions in vectors.cc + +// vectors.cc defines the vector type VEC + +template +void VectorTools::interpolate +(const Mapping&, + const DoFHandler&, + const Function&, + VEC&); +template +void VectorTools::interpolate +(const DoFHandler&, + const Function&, + VEC&); + +// Should these be instantiated for every combination of two types? +template +void VectorTools::interpolate +(const DoFHandler&, + const DoFHandler&, + const FullMatrix&, + const VEC&, + VEC&); + + +template +void VectorTools::integrate_difference +(const DoFHandler&, + const VEC&, + const Function&, + Vector&, + const Quadrature&, + const NormType&, + const Function*, + const double); +template +void VectorTools::integrate_difference +(const DoFHandler&, + const VEC&, + const Function&, + Vector&, + const Quadrature&, + const NormType&, + const Function*, + const double); +template +void VectorTools::integrate_difference +(const Mapping&, + const DoFHandler&, + const VEC&, + const Function&, + Vector&, + const Quadrature&, + const NormType&, + const Function*, + const double); +template +void VectorTools::integrate_difference +(const Mapping&, + const DoFHandler&, + const VEC&, + const Function&, + Vector&, + const Quadrature&, + const NormType&, + const Function*, + const double); + +template +void VectorTools::point_difference ( + const DoFHandler&, + const VEC&, + const Function&, + Vector&, + const Point&); + +template +double VectorTools::compute_mean_value +(const Mapping&, + const DoFHandler&, + const Quadrature&, + const VEC&, + const unsigned int); +template +double VectorTools::compute_mean_value +(const DoFHandler&, + const Quadrature&, + const VEC&, + const unsigned int); -- 2.39.5