From dcddd12cb390bd770653aa314b8de0124dba45ac Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Sat, 11 Feb 2006 00:35:33 +0000 Subject: [PATCH] Introduce the tensorial coefficient. git-svn-id: https://svn.dealii.org/trunk@12311 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-20/step-20.cc | 133 ++++++++++++++++++---------- 1 file changed, 85 insertions(+), 48 deletions(-) diff --git a/deal.II/examples/step-20/step-20.cc b/deal.II/examples/step-20/step-20.cc index 8e7c711701..53ad633f25 100644 --- a/deal.II/examples/step-20/step-20.cc +++ b/deal.II/examples/step-20/step-20.cc @@ -47,12 +47,22 @@ #include #include - // This is the only new header, - // namely the one in which the - // Raviart-Thomas finite element is - // declared: + // This is the only significant new + // header, namely the one in which + // the Raviart-Thomas finite element + // is declared: #include + // Finally, as a bonus in this + // program, we will use a tensorial + // coefficient. Since it may have a + // spatial dependence, we consider it + // a tensor-valued function. The + // following include file provides + // the ``TensorFunction'' class that + // offers such functionality: +#include + // @sect3{The ``MixedLaplaceProblem'' class template} @@ -114,7 +124,7 @@ class MixedLaplaceProblem }; - // @sect3{Right hand side, coefficient, and exact solution} + // @sect3{Right hand side, boundary values, and exact solution} // Our next task is to define the // right hand side of our problem @@ -215,6 +225,38 @@ ExactSolution::vector_value (const Point &p, + // @sect3{The permability tensor K} + +template +class Coefficient : public TensorFunction<2,dim> +{ + public: + virtual void value_list (const std::vector > &points, + std::vector > &values) const; +}; + + +template +void +Coefficient::value_list (const std::vector > &points, + std::vector > &values) const +{ + Assert (points.size() == values.size(), + ExcDimensionMismatch (points.size(), values.size())); + + for (unsigned int p=0; p MixedLaplaceProblem::MixedLaplaceProblem (const unsigned int degree) : @@ -275,61 +317,47 @@ void MixedLaplaceProblem::make_grid_and_dofs () -Tensor<1,2> extract_u (const FEValuesBase<2> &fe_values, - const unsigned int j, - const unsigned int q) +template +Tensor<1,dim> +extract_u (const FEValuesBase &fe_values, + const unsigned int i, + const unsigned int q) { - Tensor<1,2> tmp; - tmp[0] = fe_values.shape_value_component (j,q,0); - tmp[1] = fe_values.shape_value_component (j,q,1); - return tmp; -} - + Tensor<1,dim> tmp; + for (unsigned int d=0; d extract_u (const FEValuesBase<3> &fe_values, - const unsigned int j, - const unsigned int q) -{ - Tensor<1,3> tmp; - tmp[0] = fe_values.shape_value_component (j,q,0); - tmp[1] = fe_values.shape_value_component (j,q,1); - tmp[2] = fe_values.shape_value_component (j,q,2); return tmp; } - - -double extract_div_u (const FEValuesBase<2> &fe_values, - const unsigned int j, - const unsigned int q) +template +double +extract_div_u (const FEValuesBase &fe_values, + const unsigned int i, + const unsigned int q) { - return fe_values.shape_grad_component (j,q,0)[0] + - fe_values.shape_grad_component (j,q,1)[1]; -} - + double divergence = 0; + for (unsigned int d=0; d &fe_values, - const unsigned int j, - const unsigned int q) -{ - return fe_values.shape_grad_component (j,q,0)[0] + - fe_values.shape_grad_component (j,q,1)[1] + - fe_values.shape_grad_component (j,q,2)[2]; + return divergence; } + template double extract_p (const FEValuesBase &fe_values, - const unsigned int j, + const unsigned int i, const unsigned int q) { - return fe_values.shape_value_component (j,q,dim); + return fe_values.shape_value_component (i,q,dim); } + template void MixedLaplaceProblem::assemble_system () { @@ -351,14 +379,17 @@ void MixedLaplaceProblem::assemble_system () Vector local_rhs (dofs_per_cell); - const RightHandSide right_hand_side; + const RightHandSide right_hand_side; const PressureBoundaryValues pressure_boundary_values; - + const Coefficient coefficient; + std::vector rhs_values (n_q_points); std::vector boundary_values (n_face_q_points); + std::vector > Kinverse (n_q_points); std::vector local_dof_indices (dofs_per_cell); + typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); @@ -370,6 +401,8 @@ void MixedLaplaceProblem::assemble_system () right_hand_side.value_list (fe_values.get_quadrature_points(), rhs_values); + coefficient.value_list (fe_values.get_quadrature_points(), + Kinverse); for (unsigned int q=0; q::assemble_system () const double div_phi_j_u = extract_div_u (fe_values, j, q); const double phi_j_p = extract_p (fe_values, j, q); - local_matrix(i,j) += (phi_i_u * phi_j_u + local_matrix(i,j) += (phi_i_u * Kinverse[q] * phi_j_u - div_phi_i_u * phi_j_p - phi_i_p * div_phi_j_u) * fe_values.JxW(q); @@ -395,18 +428,22 @@ void MixedLaplaceProblem::assemble_system () fe_values.JxW(q); } - for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) + for (unsigned int face_no=0; + face_no::faces_per_cell; + ++face_no) if (cell->at_boundary(face_no)) { fe_face_values.reinit (cell, face_no); - pressure_boundary_values.value_list (fe_face_values.get_quadrature_points(), - boundary_values); + pressure_boundary_values + .value_list (fe_face_values.get_quadrature_points(), + boundary_values); for (unsigned int q=0; q phi_i_u = extract_u (fe_face_values, i, q); + const Tensor<1,dim> + phi_i_u = extract_u (fe_face_values, i, q); local_rhs(i) += -(phi_i_u * fe_face_values.normal_vector(q) * @@ -609,7 +646,7 @@ int main () { deallog.depth_console (0); - MixedLaplaceProblem<2> mixed_laplace_problem (1); + MixedLaplaceProblem<2> mixed_laplace_problem (0); mixed_laplace_problem.run (); return 0; -- 2.39.5