From 6d95c72baafdaa0d501e60ef2ef173aa3a1c41e0 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Sun, 28 Feb 2016 11:38:21 +0100 Subject: [PATCH] extend MatrixTools::create_boundary_mass_matrix() for complex algebra Most of the modifications are change of the multiplication order from double * std::complex to std::complex * double. Instantiate create_mass_matrix() and create_boundary_mass_matrix() for complex case. --- source/numerics/matrix_tools.cc | 103 +++++++++---------- source/numerics/matrix_tools.inst.in | 141 +++++++++++++++++++++++++++ 2 files changed, 194 insertions(+), 50 deletions(-) diff --git a/source/numerics/matrix_tools.cc b/source/numerics/matrix_tools.cc index abeaaae37a..44f2e1dbc1 100644 --- a/source/numerics/matrix_tools.cc +++ b/source/numerics/matrix_tools.cc @@ -617,7 +617,7 @@ namespace MatrixCreator Scratch() {} }; - template + template struct CopyData { CopyData() {}; @@ -628,12 +628,12 @@ namespace MatrixCreator std::vector dofs; std::vector > dof_is_on_face; typename DoFHandlerType::active_cell_iterator cell; - std::vector > cell_matrix; - std::vector > cell_vector; + std::vector > cell_matrix; + std::vector > cell_vector; }; - template - CopyData::CopyData(CopyData const &data) : + template + CopyData::CopyData(CopyData const &data) : dofs_per_cell(data.dofs_per_cell), dofs(data.dofs), dof_is_on_face(data.dof_is_on_face), @@ -876,7 +876,7 @@ namespace MatrixCreator create_boundary_mass_matrix_1 (typename DoFHandler::active_cell_iterator const &cell, MatrixCreator::internal::AssemblerBoundary::Scratch const &, MatrixCreator::internal::AssemblerBoundary::CopyData > ©_data, + spacedim>,number > ©_data, Mapping const &mapping, FiniteElement const &fe, Quadrature const &q, @@ -935,9 +935,9 @@ namespace MatrixCreator if (boundary_functions.find(cell->face(face)->boundary_id()) != boundary_functions.end()) { - copy_data.cell_matrix.push_back(FullMatrix (copy_data.dofs_per_cell, + copy_data.cell_matrix.push_back(FullMatrix (copy_data.dofs_per_cell, copy_data.dofs_per_cell)); - copy_data.cell_vector.push_back(Vector (copy_data.dofs_per_cell)); + copy_data.cell_vector.push_back(Vector (copy_data.dofs_per_cell)); fe_values.reinit (cell, face); if (fe_is_system) @@ -1005,14 +1005,14 @@ namespace MatrixCreator == fe.system_to_component_index(i).first) { copy_data.cell_matrix.back()(i,j) - += weight + += coefficient_vector_values[point](fe.system_to_component_index(i).first) + * weight * fe_values.shape_value(j,point) - * fe_values.shape_value(i,point) - * coefficient_vector_values[point](fe.system_to_component_index(i).first); + * fe_values.shape_value(i,point); } } - copy_data.cell_vector.back()(i) += fe_values.shape_value(i,point) - * rhs_values_system[point](component_mapping[fe.system_to_component_index(i).first]) + copy_data.cell_vector.back()(i) += rhs_values_system[point](component_mapping[fe.system_to_component_index(i).first]) + * fe_values.shape_value(i,point) * weight; } else @@ -1021,12 +1021,13 @@ namespace MatrixCreator { for (unsigned int j=0; j void copy_boundary_mass_matrix_1(MatrixCreator::internal::AssemblerBoundary::CopyData > const ©_data, + spacedim>,number > const ©_data, std::map*> const &boundary_functions, std::vector const &dof_to_boundary_mapping, SparseMatrix &matrix, @@ -1145,7 +1146,7 @@ namespace MatrixCreator create_boundary_mass_matrix_1<1,3,float> (DoFHandler<1,3>::active_cell_iterator const &/*cell*/, MatrixCreator::internal::AssemblerBoundary::Scratch const &, MatrixCreator::internal::AssemblerBoundary::CopyData > &/*copy_data*/, + 3>,float > &/*copy_data*/, Mapping<1,3> const &, FiniteElement<1,3> const &, Quadrature<0> const &, @@ -1161,7 +1162,7 @@ namespace MatrixCreator create_boundary_mass_matrix_1<1,3,double> (DoFHandler<1,3>::active_cell_iterator const &/*cell*/, MatrixCreator::internal::AssemblerBoundary::Scratch const &, MatrixCreator::internal::AssemblerBoundary::CopyData > &/*copy_data*/, + 3>,double > &/*copy_data*/, Mapping<1,3> const &, FiniteElement<1,3> const &, Quadrature<0> const &, @@ -1177,7 +1178,7 @@ namespace MatrixCreator create_boundary_mass_matrix_1<1,3,long double> (DoFHandler<1,3>::active_cell_iterator const &/*cell*/, MatrixCreator::internal::AssemblerBoundary::Scratch const &, MatrixCreator::internal::AssemblerBoundary::CopyData > &/*copy_data*/, + 3>,long double > &/*copy_data*/, Mapping<1,3> const &, FiniteElement<1,3> const &, Quadrature<0> const &, @@ -1238,19 +1239,19 @@ namespace MatrixCreator AssertDimension (n_components, component_mapping.size()); MatrixCreator::internal::AssemblerBoundary::Scratch scratch; - MatrixCreator::internal::AssemblerBoundary::CopyData > copy_data; + MatrixCreator::internal::AssemblerBoundary::CopyData,number > copy_data; WorkStream::run(dof.begin_active(),dof.end(), static_cast::active_cell_iterator const &,MatrixCreator::internal::AssemblerBoundary::Scratch const &, - MatrixCreator::internal::AssemblerBoundary::CopyData > &)> > + MatrixCreator::internal::AssemblerBoundary::CopyData,number > &)> > (std_cxx11::bind(&create_boundary_mass_matrix_1,std_cxx11::_1,std_cxx11::_2, std_cxx11::_3, std_cxx11::cref(mapping),std_cxx11::cref(fe),std_cxx11::cref(q), std_cxx11::cref(boundary_functions),coefficient, std_cxx11::cref(component_mapping))), static_cast > const &)> > (std_cxx11::bind( + ::CopyData,number> const &)> > (std_cxx11::bind( ©_boundary_mass_matrix_1, std_cxx11::_1, std_cxx11::cref(boundary_functions), @@ -1272,7 +1273,7 @@ namespace MatrixCreator &cell, MatrixCreator::internal::AssemblerBoundary::Scratch const &, MatrixCreator::internal::AssemblerBoundary - ::CopyData > ©_data, + ::CopyData,number > ©_data, hp::MappingCollection const &mapping, hp::FECollection const &fe_collection, hp::QCollection const &q, @@ -1329,9 +1330,9 @@ namespace MatrixCreator const FEFaceValues &fe_values = x_fe_values.get_present_fe_values (); - copy_data.cell_matrix.push_back(FullMatrix (copy_data.dofs_per_cell, + copy_data.cell_matrix.push_back(FullMatrix (copy_data.dofs_per_cell, copy_data.dofs_per_cell)); - copy_data.cell_vector.push_back(Vector (copy_data.dofs_per_cell)); + copy_data.cell_vector.push_back(Vector (copy_data.dofs_per_cell)); if (fe_is_system) // FE has several components @@ -1361,12 +1362,13 @@ namespace MatrixCreator { const double u = fe_values.shape_value(j,point); copy_data.cell_matrix.back()(i,j) - += (u * v * weight * coefficient_values[point]); + += (coefficient_values[point] * u * v * weight); } - copy_data.cell_vector.back()(i) += v * - rhs_values_system[point]( - component_mapping[fe.system_to_component_index(i).first]) * weight; + copy_data.cell_vector.back()(i) += rhs_values_system[point]( + component_mapping[fe.system_to_component_index(i).first]) + * weight + * v; } } } @@ -1390,10 +1392,11 @@ namespace MatrixCreator { const double u = fe_values.shape_value(j,point); copy_data.cell_matrix.back()(i,j) += - (u * v * weight * coefficient_vector_values[point](component_i)); + (coefficient_vector_values[point](component_i) * u * v * weight); } - copy_data.cell_vector.back()(i) += v * - rhs_values_system[point](component_mapping[component_i]) * weight; + copy_data.cell_vector.back()(i) += rhs_values_system[point](component_mapping[component_i]) + * v + * weight; } } } @@ -1412,9 +1415,9 @@ namespace MatrixCreator const double u = fe_values.shape_value(j,point); copy_data.cell_matrix.back()(i,j) += (u * v * weight); } - copy_data.cell_vector.back()(i) += v * - rhs_values_system[point]( + copy_data.cell_vector.back()(i) += rhs_values_system[point]( fe.system_to_component_index(i).first) * + v * weight; } } @@ -1440,10 +1443,10 @@ namespace MatrixCreator for (unsigned int j=0; j void copy_hp_boundary_mass_matrix_1(MatrixCreator::internal::AssemblerBoundary - ::CopyData > const ©_data, + ::CopyData,number > const ©_data, std::map*> const &boundary_functions, std::vector const &dof_to_boundary_mapping, SparseMatrix &matrix, @@ -1547,8 +1550,8 @@ namespace MatrixCreator double max_diag_entry = 0; for (unsigned int i=0; i max_diag_entry) - max_diag_entry = std::fabs(copy_data.cell_matrix[pos](i,i)); + if (std::abs(copy_data.cell_matrix[pos](i,i)) > max_diag_entry) + max_diag_entry = std::abs(copy_data.cell_matrix[pos](i,i)); #endif for (unsigned int i=0; i > copy_data; + MatrixCreator::internal::AssemblerBoundary::CopyData,number > copy_data; WorkStream::run(dof.begin_active(),dof.end(), static_cast::active_cell_iterator const &,MatrixCreator::internal::AssemblerBoundary::Scratch const &, - MatrixCreator::internal::AssemblerBoundary::CopyData > &)> > + MatrixCreator::internal::AssemblerBoundary::CopyData,number > &)> > (std_cxx11::bind( &create_hp_boundary_mass_matrix_1,std_cxx11::_1,std_cxx11::_2, std_cxx11::_3, std_cxx11::cref(mapping),std_cxx11::cref(fe_collection),std_cxx11::cref(q), std_cxx11::cref(boundary_functions),coefficient, std_cxx11::cref(component_mapping))), static_cast > const &)> > ( + ::CopyData,number > const &)> > ( std_cxx11::bind( ©_hp_boundary_mass_matrix_1, std_cxx11::_1, std_cxx11::cref(boundary_functions), diff --git a/source/numerics/matrix_tools.inst.in b/source/numerics/matrix_tools.inst.in index 59073c0117..6ad857c2a4 100644 --- a/source/numerics/matrix_tools.inst.in +++ b/source/numerics/matrix_tools.inst.in @@ -93,6 +93,47 @@ for (scalar: REAL_SCALARS; deal_II_dimension : DIMENSIONS; deal_II_space_dimensi #endif } +for (scalar: COMPLEX_SCALARS; deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) + { +#if deal_II_dimension <= deal_II_space_dimension + // non-hp version of create_mass_matrix + template + void MatrixCreator::create_mass_matrix + (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function * const coefficient, + const ConstraintMatrix &constraints); + template + void MatrixCreator::create_mass_matrix + (const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function * const coefficient, + const ConstraintMatrix &constraints); + template + void MatrixCreator::create_mass_matrix + (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function * const coefficient, + const ConstraintMatrix &constraints); + template + void MatrixCreator::create_mass_matrix + (const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function * const coefficient, + const ConstraintMatrix &constraints); +#endif + } + for (scalar: REAL_SCALARS; deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) { @@ -146,6 +187,57 @@ for (scalar: REAL_SCALARS; deal_II_dimension : DIMENSIONS; deal_II_space_dimensi #endif } +for (scalar: COMPLEX_SCALARS; deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) + { +#if deal_II_dimension <= deal_II_space_dimension + template + void MatrixCreator::create_boundary_mass_matrix + (const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const std::map*> &rhs, + Vector &rhs_vector, + std::vector &dof_to_boundary_mapping, + const Function * const a, + std::vector); + + template + void MatrixCreator::create_boundary_mass_matrix + (const Mapping &, + const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const std::map*> &rhs, + Vector &rhs_vector, + std::vector &dof_to_boundary_mapping, + const Function * const a, + std::vector); + + template + void + MatrixCreator::create_boundary_mass_matrix + (const hp::MappingCollection&, + const hp::DoFHandler&, + const hp::QCollection&, + SparseMatrix&, + const std::map*>&, + Vector&, + std::vector&, + const Function * const, + std::vector); + + template + void MatrixCreator::create_boundary_mass_matrix + (const hp::DoFHandler&, + const hp::QCollection&, + SparseMatrix&, + const std::map*>&, + Vector&, + std::vector&, + const Function * const, + std::vector); +#endif + } //TODO[SP]: replace by // where applicable and move to codimension cases above also when applicable @@ -199,6 +291,55 @@ for (scalar : REAL_SCALARS; deal_II_dimension : DIMENSIONS; deal_II_space_dimens #endif } +for (scalar : COMPLEX_SCALARS; deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) + { +// hp versions of functions +#if deal_II_dimension <= deal_II_space_dimension + template + void MatrixCreator::create_mass_matrix + (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const Function * const coefficient, + const ConstraintMatrix &constraints); + + template + void MatrixCreator::create_mass_matrix + (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function * const coefficient, + const ConstraintMatrix &constraints); + +#endif + +#if deal_II_dimension == deal_II_space_dimension + + + template + void MatrixCreator::create_mass_matrix + (const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const Function * const coefficient, + const ConstraintMatrix &constraints); + + template + void MatrixCreator::create_mass_matrix + (const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function * const coefficient, + const ConstraintMatrix &constraints); + +#endif + } for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) { -- 2.39.5