From 06b2946e1641d445160f1036f77dfa2d9063d9c7 Mon Sep 17 00:00:00 2001 From: guido Date: Wed, 24 Feb 1999 13:19:48 +0000 Subject: [PATCH] Preparation for systems and mixed discretizations git-svn-id: https://svn.dealii.org/trunk@893 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/Makefile | 2 +- deal.II/deal.II/include/fe/fe.h | 10 +- deal.II/deal.II/include/fe/fe_system.h | 24 +++ deal.II/deal.II/include/numerics/vectors.h | 16 ++ deal.II/deal.II/source/fe/fe.cc | 9 +- .../deal.II/source/fe/fe_lib.criss_cross.cc | 5 +- .../deal.II/source/fe/fe_linear_mapping.cc | 11 +- deal.II/deal.II/source/fe/fe_system.cc | 130 ++++++++---- deal.II/deal.II/source/fe/fe_values.cc | 33 ++- deal.II/deal.II/source/numerics/data_io.cc | 21 +- deal.II/deal.II/source/numerics/vectors.cc | 188 ++++++++++++++++++ 11 files changed, 383 insertions(+), 66 deletions(-) diff --git a/deal.II/deal.II/Makefile b/deal.II/deal.II/Makefile index e277543fe0..fa58224016 100644 --- a/deal.II/deal.II/Makefile +++ b/deal.II/deal.II/Makefile @@ -23,7 +23,7 @@ examples-clean: source-clean: cd source ; make clean TAGS: - etags include/*/* source/*/* + etags include/*/*.h source/*/*.cc .PHONY: all deal.II examples TAGS diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h index 7d29d97057..b6d47ef7f4 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -92,8 +92,9 @@ class FiniteElementData */ const unsigned int n_transform_functions; + /** - * Number of components. + * Number of components and dimension of the image space. */ const unsigned int n_components; @@ -853,7 +854,12 @@ class FiniteElement : public FiniteElementBase { * * The function assumes that the fields * already have the right number of - * elements. + * elements. It has to be + * guaranteed, that fields that are + * not requested for update are not changed. + * This also means, that these + * fields have to be filled with + * the correct values beforehand. * * This function is more or less an * interface to the #FEValues# class and diff --git a/deal.II/deal.II/include/fe/fe_system.h b/deal.II/deal.II/include/fe/fe_system.h index d702635683..00c152ec23 100644 --- a/deal.II/deal.II/include/fe/fe_system.h +++ b/deal.II/deal.II/include/fe/fe_system.h @@ -316,6 +316,23 @@ class FESystem : public FiniteElement const vector > &unit_points, vector > &normal_vectors) const; + /** + * Implementation of the corresponding function of #FiniteElement#. + */ + virtual void fill_fe_values (const DoFHandler::cell_iterator &cell, + const vector > &unit_points, + vector > &jacobians, + const bool compute_jacobians, + vector > &jacobians_grad, + const bool compute_jacobians_grad, + vector > &support_points, + const bool compute_support_points, + vector > &q_points, + const bool compute_q_points, + const dFMatrix &shape_values_transform, + const vector > > &shape_grad_transform, + const Boundary &boundary) const; + /** * Number of different base * elements of this object. @@ -425,6 +442,10 @@ class FESystem : public FiniteElement * the #.h# file. */ void initialize(); + /** + *Exception. + */ + DeclException0(ExcElementTransformNotEqual); }; @@ -478,6 +499,9 @@ FESystem::FESystem (const FE1 &fe1, const unsigned int n1, base_elements(2), component_to_base_table(n_components) { + Assert(fe1.n_transform_functions == fe2.n_transform_functions, + ExcElementTransformNotEqual()); + base_elements[0] = ElementPair(new FE1, n1); base_elements[0].first -> subscribe (); base_elements[1] = ElementPair(new FE2, n2); diff --git a/deal.II/deal.II/include/numerics/vectors.h b/deal.II/deal.II/include/numerics/vectors.h index 8a32f053c3..90720acd88 100644 --- a/deal.II/deal.II/include/numerics/vectors.h +++ b/deal.II/deal.II/include/numerics/vectors.h @@ -19,6 +19,7 @@ template class Boundary; template class StraightBoundary; class ConstraintMatrix; class dVector; +class VectorFunction; /** @@ -376,6 +377,7 @@ class VectorTools { map &boundary_values); /** + * Compute the error of the finite element solution. * Integrate the difference between * a finite element function and * the reference function, which @@ -393,6 +395,20 @@ class VectorTools { const NormType &norm, const Boundary &boundary=StraightBoundary()); + /** + * Compute the error for the solution of a system. + * See the other #integrate_difference#. + */ + static void integrate_difference (const DoFHandler &dof, + const dVector &fe_function, + const VectorFunction &exact_solution, + dVector &difference, + const Quadrature &q, + const FiniteElement &fe, + const NormType &norm, + const Boundary &boundary); + + /** * Exception */ diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index bf65e97833..32c438a615 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -99,7 +99,6 @@ FiniteElementData::FiniteElementData (const unsigned int dofs_per_vertex, }; - template bool FiniteElementData::operator== (const FiniteElementData &f) const { @@ -241,11 +240,13 @@ void FiniteElement<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell, const dFMatrix &, const vector > > &, const Boundary<1> &boundary) const { - Assert (jacobians.size() == unit_points.size(), + Assert ((!compute_jacobians) || (jacobians.size() == unit_points.size()), ExcWrongFieldDimension(jacobians.size(), unit_points.size())); - Assert (q_points.size() == unit_points.size(), + Assert ((!compute_jacobians_grad) || (jacobians_grad.size() == unit_points.size()), + ExcWrongFieldDimension(jacobians_grad.size(), unit_points.size())); + Assert ((!compute_q_points) || (q_points.size() == unit_points.size()), ExcWrongFieldDimension(q_points.size(), unit_points.size())); - Assert (support_points.size() == total_dofs, + Assert ((!compute_support_points) || (support_points.size() == total_dofs), ExcWrongFieldDimension(support_points.size(), total_dofs)); diff --git a/deal.II/deal.II/source/fe/fe_lib.criss_cross.cc b/deal.II/deal.II/source/fe/fe_lib.criss_cross.cc index 40f21263bd..6000e8c3c3 100644 --- a/deal.II/deal.II/source/fe/fe_lib.criss_cross.cc +++ b/deal.II/deal.II/source/fe/fe_lib.criss_cross.cc @@ -939,7 +939,7 @@ void FECrissCross::fill_fe_values (const DoFHandler::cell_iterator &ce vector > &jacobians_grad, const bool compute_jacobians_grad, vector > &support_points, - const bool, + const bool compute_support_points, vector > &q_points, const bool compute_q_points, const dFMatrix &shape_values_transform, @@ -957,7 +957,8 @@ void FECrissCross::fill_fe_values (const DoFHandler::cell_iterator &ce // we need the support points in any // way, wanted or not by the user - get_support_points (cell, boundary, support_points); + if (compute_support_points) + get_support_points (cell, boundary, support_points); if (compute_q_points) { diff --git a/deal.II/deal.II/source/fe/fe_linear_mapping.cc b/deal.II/deal.II/source/fe/fe_linear_mapping.cc index 0d0acd4397..bef427b76a 100644 --- a/deal.II/deal.II/source/fe/fe_linear_mapping.cc +++ b/deal.II/deal.II/source/fe/fe_linear_mapping.cc @@ -567,12 +567,15 @@ void FELinearMapping::fill_fe_values (const DoFHandler::cell_iterator const bool compute_q_points, const dFMatrix &shape_values_transform, const vector > > &/*shape_grad_transform*/, - const Boundary &boundary) const { - Assert (jacobians.size() == unit_points.size(), + const Boundary &boundary) const +{ + Assert ((!compute_jacobians) || (jacobians.size() == unit_points.size()), ExcWrongFieldDimension(jacobians.size(), unit_points.size())); - Assert (q_points.size() == unit_points.size(), + Assert ((!compute_jacobians_grad) || (jacobians_grad.size() == unit_points.size()), + ExcWrongFieldDimension(jacobians_grad.size(), unit_points.size())); + Assert ((!compute_q_points) || (q_points.size() == unit_points.size()), ExcWrongFieldDimension(q_points.size(), unit_points.size())); - Assert (support_points.size() == total_dofs, + Assert ((!compute_support_points) || (support_points.size() == total_dofs), ExcWrongFieldDimension(support_points.size(), total_dofs)); diff --git a/deal.II/deal.II/source/fe/fe_system.cc b/deal.II/deal.II/source/fe/fe_system.cc index 94991e57d4..6ea9c6f79a 100644 --- a/deal.II/deal.II/source/fe/fe_system.cc +++ b/deal.II/deal.II/source/fe/fe_system.cc @@ -1,5 +1,5 @@ /* $Id$ */ -/* Copyright W. Bangerth, University of Heidelberg, 1990 */ +/* Copyright W. Bangerth, G. Kanschat University of Heidelberg, 1990 */ #include @@ -196,7 +196,7 @@ FESystem<1>::multiply_dof_numbers (const FiniteElementData<1> &fe1, { return FiniteElementData<1> (fe1.dofs_per_vertex * N1 + fe2.dofs_per_vertex * N2 , fe1.dofs_per_line * N1 + fe2.dofs_per_line * N2 , - fe1.n_transform_functions + fe2.n_transform_functions , + fe1.n_transform_functions, fe1.n_components * N1 + fe2.n_components * N2 ); }; @@ -227,7 +227,7 @@ FESystem<2>::multiply_dof_numbers (const FiniteElementData<2> &fe1, return FiniteElementData<2> (fe1.dofs_per_vertex * N1 + fe2.dofs_per_vertex * N2 , fe1.dofs_per_line * N1 + fe2.dofs_per_line * N2 , fe1.dofs_per_quad * N1 + fe2.dofs_per_quad * N2 , - fe1.n_transform_functions + fe2.n_transform_functions , + fe1.n_transform_functions, fe1.n_components * N1 + fe2.n_components * N2 ); }; @@ -374,82 +374,126 @@ void FESystem::get_local_mass_matrix (const DoFHandler::cell_iterator template -double FESystem::shape_value_transform (const unsigned int /*i*/, - const Point &/*p*/) const +double FESystem::shape_value_transform (const unsigned int i, + const Point &p) const { - Assert(false, ExcNotImplemented()); - return 0.; + return base_elements[0].first->shape_value_transform(i,p); }; template -Tensor<1,dim> FESystem::shape_grad_transform (const unsigned int /*i*/, - const Point &/*p*/) const +Tensor<1,dim> FESystem::shape_grad_transform (const unsigned int i, + const Point &p) const { - Assert(false, ExcNotImplemented()); - return Tensor<1,dim>(); - -// return base_element->shape_grad_transform (i, p); + return base_elements[0].first->shape_grad_transform (i, p); }; template -void FESystem::get_face_jacobians (const DoFHandler::face_iterator &/*face*/, - const Boundary &/*boundary*/, - const vector > &/*unit_points*/, - vector &/*face_jacobi_determinants*/) const +void FESystem::get_face_jacobians (const DoFHandler::face_iterator &face, + const Boundary &boundary, + const vector > &unit_points, + vector &face_jacobi_determinants) const { - Assert(false, ExcNotImplemented()); - -// base_element->get_face_jacobians (face, boundary, unit_points, face_jacobi_determinants); + base_elements[0].first->get_face_jacobians (face, boundary, unit_points, + face_jacobi_determinants); }; template -void FESystem::get_subface_jacobians (const DoFHandler::face_iterator &/*face*/, - const unsigned int /*subface_no*/, - const vector > &/*unit_points*/, - vector &/*face_jacobi_determinants*/) const +void FESystem::get_subface_jacobians (const DoFHandler::face_iterator &face, + const unsigned int subface_no, + const vector > &unit_points, + vector &face_jacobi_determinants) const { - Assert(false, ExcNotImplemented()); - -// base_element->get_subface_jacobians (face, subface_no, unit_points, face_jacobi_determinants); + base_elements[0].first->get_subface_jacobians (face, subface_no, unit_points, + face_jacobi_determinants); }; template -void FESystem::get_normal_vectors (const DoFHandler::cell_iterator &/*cell*/, - const unsigned int /*face_no*/, - const Boundary &/*boundary*/, - const vector > &/*unit_points*/, - vector > &/*normal_vectors*/) const +void FESystem::get_normal_vectors (const DoFHandler::cell_iterator &cell, + const unsigned int face_no, + const Boundary &boundary, + const vector > &unit_points, + vector > &normal_vectors) const { - Assert(false, ExcNotImplemented()); - -// base_element->get_normal_vectors (cell, face_no, boundary, unit_points, normal_vectors); + base_elements[0].first->get_normal_vectors (cell, face_no, boundary, unit_points, + normal_vectors); }; template -void FESystem::get_normal_vectors (const DoFHandler::cell_iterator &/*cell*/, - const unsigned int /*face_no*/, - const unsigned int /*subface_no*/, - const vector > &/*unit_points*/, - vector > &/*normal_vectors*/) const +void FESystem::get_normal_vectors (const DoFHandler::cell_iterator &cell, + const unsigned int face_no, + const unsigned int subface_no, + const vector > &unit_points, + vector > &normal_vectors) const { - Assert(false, ExcNotImplemented()); - -// base_element->get_normal_vectors (cell, face_no, subface_no, unit_points, normal_vectors); + base_elements[0].first->get_normal_vectors (cell, face_no, subface_no, unit_points, + normal_vectors); }; +template +void +FESystem::fill_fe_values (const DoFHandler::cell_iterator &cell, + const vector > &unit_points, + vector > &jacobians, + const bool compute_jacobians, + vector > &jacobians_grad, + const bool compute_jacobians_grad, + vector > &support_points, + const bool compute_support_points, + vector > &q_points, + const bool compute_q_points, + const dFMatrix &shape_values_transform, + const vector > > &shape_grad_transform, + const Boundary &boundary) const +{ + vector > supp(base_elements[0].first->total_dofs); + base_elements[0].first->fill_fe_values (cell, unit_points, jacobians, compute_jacobians, + jacobians_grad, compute_jacobians_grad, + support_points, compute_support_points, + q_points, compute_q_points, + shape_values_transform, shape_grad_transform, boundary); + + if (compute_support_points) + { + unsigned component = 0; + + for (unsigned m=0 ; m < element_multiplicity(0) ; ++ m) + { + for (unsigned i=0 ; i < base_element(0).total_dofs ; ++i) + support_points[component_to_system_index(component,i)] = supp[i]; + ++component; + } + for (unsigned base=1 ; base < n_base_elements() ; ++base) + { + supp.resize(base_elements[base].first->total_dofs); + base_elements[base].first->fill_fe_values (cell, unit_points, jacobians, false, + jacobians_grad, false, + supp, true, + q_points, false, + shape_values_transform, shape_grad_transform, boundary); + + for (unsigned m=0 ; m < element_multiplicity(base) ; ++ m) + { + for (unsigned i=0 ; i < base_element(base).total_dofs ; ++i) + support_points[component_to_system_index(component,i)] = supp[i]; + ++component; + } + } + } +} + // explicit instantiations template class FESystem; diff --git a/deal.II/deal.II/source/fe/fe_values.cc b/deal.II/deal.II/source/fe/fe_values.cc index 8b51c2e287..870594cd7a 100644 --- a/deal.II/deal.II/source/fe/fe_values.cc +++ b/deal.II/deal.II/source/fe/fe_values.cc @@ -88,11 +88,42 @@ void FEValuesBase::get_function_values (const dVector &fe_function, }; +template +void FEValuesBase::get_function_values (const dVector &fe_function, + vector< vector > &values) const +{ + Assert (fe->n_components == values.size(), + ExcWrongNoOfComponents()); + Assert (selected_datasetget_dof_values (fe_function, dof_values); + + // initialize with zero + for (unsigned i=0;isystem_to_component_index(shape_func).first][point] + += (dof_values(shape_func) * shape_values[selected_dataset](shape_func, point)); +}; + + template const Tensor<1,dim> & FEValuesBase::shape_grad (const unsigned int i, - const unsigned int j) const { + const unsigned int j) const +{ Assert (i::write_gnuplot (ostream &out, unsigned int accuracy) const FEValues fe(dofs->get_fe(), points, UpdateFlags(update_q_points)); const StraightBoundary boundary; - vector *values = new vector [data.size()]; + vector< vector > > + values (data.size(), vector< vector >(dofs->get_fe().n_components, vector(points.n_quadrature_points))); for (cell=dofs->begin_active(); cell!=endc; ++cell) { @@ -551,7 +552,7 @@ void DataOut::write_gnuplot (ostream &out, unsigned int accuracy) const for (unsigned i=0; i::write_gnuplot (ostream &out, unsigned int accuracy) const Point pt = fe.quadrature_point(supp_pt); out << pt << " "; for (unsigned int i=0; i!=data.size(); ++i) - out << values[i][supp_pt] - << ' '; + for (unsigned int j=0; j < dofs->get_fe().n_components; ++j) + out << values[i][j][supp_pt] + << ' '; out << endl; }; @@ -589,8 +591,9 @@ void DataOut::write_gnuplot (ostream &out, unsigned int accuracy) const out << pt << " "; for (unsigned int i=0; i!=data.size(); ++i) - out << values[i][supp_pt] - << ' '; + for (unsigned int j=0; j < dofs->get_fe().n_components; ++j) + out << values[i][j][supp_pt] + << ' '; out << endl; } out << endl; @@ -621,8 +624,9 @@ void DataOut::write_gnuplot (ostream &out, unsigned int accuracy) const out << pt << " "; for (unsigned int i=0; i!=data.size(); ++i) - out << values[i][supp_pt] - << ' '; + for (unsigned int j=0; j < dofs->get_fe().n_components; ++j) + out << values[i][j][supp_pt] + << ' '; out << endl; } out << endl; @@ -638,7 +642,6 @@ void DataOut::write_gnuplot (ostream &out, unsigned int accuracy) const } out << endl; } - delete [] values; AssertThrow (out, ExcIO()); } diff --git a/deal.II/deal.II/source/numerics/vectors.cc b/deal.II/deal.II/source/numerics/vectors.cc index 9d8213c380..58fe08915e 100644 --- a/deal.II/deal.II/source/numerics/vectors.cc +++ b/deal.II/deal.II/source/numerics/vectors.cc @@ -555,5 +555,193 @@ void VectorTools::integrate_difference (const DoFHandler &dof, }; }; +template +void VectorTools::integrate_difference (const DoFHandler &dof, + const dVector &fe_function, + const TensorFunction<1,dim> &exact_solution, + dVector &difference, + const Quadrature &q, + const FiniteElement &fe, + const NormType &norm, + const Boundary &boundary) { + Assert (fe == dof.get_fe(), ExcInvalidFE()); + + difference.reinit (dof.get_tria().n_active_cells()); + + UpdateFlags update_flags = UpdateFlags (update_q_points | + update_JxW_values); + if ((norm==H1_seminorm) || (norm==H1_norm)) + update_flags = UpdateFlags (update_flags | update_gradients); + FEValues fe_values(fe, q, update_flags); + + // loop over all cells + 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, boundary); + + switch (norm) + { + case mean: + case L1_norm: + case L2_norm: + case Linfty_norm: + case H1_norm: + { + // we need the finite element + // function \psi at the different + // integration points. Compute + // it like this: + // \psi(x_j)=\sum_i v_i \phi_i(x_j) + // with v_i the nodal values of the + // fe_function and \phi_i(x_j) the + // matrix of the trial function + // values at the integration point + // x_j. Then the vector + // of the \psi(x_j) is v*Phi with + // v being the vector of nodal + // values on this cell and Phi + // the matrix. + // + // we then need the difference: + // reference_function(x_j)-\psi_j + // and assign that to the vector + // \psi. + const unsigned int n_q_points = q.n_quadrature_points; + vector psi (n_q_points); + + // in praxi: first compute + // exact fe_function vector + exact_solution.value_list (fe_values.get_quadrature_points(), + psi); + // then subtract finite element + // fe_function + if (true) + { + vector< vector > function_values (fe.n_components, vector(n_q_points, 0)); + fe_values.get_function_values (fe_function, function_values); + + transform (psi.begin(), psi.end(), + function_values.begin(), + psi.begin(), + minus()); + }; + + // for L1_norm and Linfty_norm: + // take absolute + // value, for the L2_norm take + // square of psi + switch (norm) + { + case mean: + break; + case L1_norm: + case Linfty_norm: + transform (psi.begin(), psi.end(), + psi.begin(), ptr_fun(fabs)); + break; + case L2_norm: + case H1_norm: + transform (psi.begin(), psi.end(), + psi.begin(), ptr_fun(sqr)); + break; + default: + Assert (false, ExcNotImplemented()); + }; + + // ok, now we have the integrand, + // let's compute the integral, + // which is + // sum_j psi_j JxW_j + // (or |psi_j| or |psi_j|^2 + switch (norm) + { + case mean: + case L1_norm: + diff = inner_product (psi.begin(), psi.end(), + fe_values.get_JxW_values().begin(), + 0.0); + break; + case L2_norm: + case H1_norm: + diff = sqrt(inner_product (psi.begin(), psi.end(), + fe_values.get_JxW_values().begin(), + 0.0)); + break; + case Linfty_norm: + diff = *max_element (psi.begin(), psi.end()); + break; + default: + Assert (false, ExcNotImplemented()); + }; + + // note: the H1_norm uses the result + // of the L2_norm and control goes + // over to the next case statement! + if (norm != H1_norm) + break; + }; + + case H1_seminorm: + { + // note: the computation of the + // H1_norm starts at the previous + // case statement, but continues + // here! + + // for H1_norm: re-square L2_norm. + diff = sqr(diff); + + // same procedure as above, but now + // psi is a vector of gradients + const unsigned int n_q_points = q.n_quadrature_points; + vector > psi (n_q_points); + + // in praxi: first compute + // exact fe_function vector + exact_solution.gradient_list (fe_values.get_quadrature_points(), + psi); + + // then subtract finite element + // fe_function + if (true) + { + vector > function_grads (n_q_points, Tensor<1,dim>()); + fe_values.get_function_grads (fe_function, function_grads); + + transform (psi.begin(), psi.end(), + function_grads.begin(), + psi.begin(), + minus >()); + }; + // take square of integrand + vector psi_square (psi.size(), 0.0); + for (unsigned int i=0; i; -- 2.39.5