From: Martin Kronbichler Date: Fri, 7 Jul 2017 20:19:11 +0000 (+0200) Subject: Use tensor product matrix in MappingQGeneric. X-Git-Tag: v9.0.0-rc1~1419^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=93d51c654e9d6b92a4157a7c39ad24d7dcd64b17;p=dealii.git Use tensor product matrix in MappingQGeneric. --- diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc index 292af5d8b9..fe1c656089 100644 --- a/source/fe/mapping_q_generic.cc +++ b/source/fe/mapping_q_generic.cc @@ -22,6 +22,7 @@ #include #include #include +#include #include #include #include @@ -32,6 +33,7 @@ #include #include #include +#include #include #include @@ -852,9 +854,6 @@ namespace internal dealii::Table<2,double> compute_laplace_vector(const unsigned int polynomial_degree) { - dealii::Table<2,double> lvs; - - Assert(lvs.n_rows()==0, ExcInternalError()); Assert(dim==2 || dim==3, ExcNotImplemented()); // for degree==1, we shouldn't have to compute any support points, since all @@ -867,8 +866,10 @@ namespace internal 4+4*(polynomial_degree-1) : 8+12*(polynomial_degree-1)+6*(polynomial_degree-1)*(polynomial_degree-1)); + dealii::Table<2,double> lvs(n_inner, n_outer); - // compute the shape gradients at the quadrature points on the unit cell + // compute the shape gradients at the quadrature points on the unit + // cell const QGauss quadrature(polynomial_degree+1); const unsigned int n_q_points=quadrature.size(); @@ -877,6 +878,10 @@ namespace internal n_q_points); quadrature_data.compute_shape_function_values(quadrature.get_points()); +#ifndef DEAL_II_WITH_LAPACK + + // Slow implementation: matrix-based method + // Compute the stiffness matrix of the inner dofs FullMatrix S(n_inner); for (unsigned int point=0; point T(n_inner, n_outer); for (unsigned int point=0; point gauss_1d(polynomial_degree+1); + FE_Q<1> fe_1d(polynomial_degree); + AlignedVector value_1d((polynomial_degree-1)*(polynomial_degree+1)); + AlignedVector derivative_1d((polynomial_degree-1)*(polynomial_degree+1)); + for (unsigned int i=0; i mass_1d(polynomial_degree-1, polynomial_degree-1); + FullMatrix lapl_1d(mass_1d); + for (unsigned int i=0; i tensor_product_matrix; + tensor_product_matrix.reinit(mass_1d, lapl_1d); + + internal::EvaluatorTensorProduct + eval_rhs(value_1d, derivative_1d, value_1d, + polynomial_degree-2, polynomial_degree+1); + + std::vector tmp_rhs(dim*n_q_points), rl1(n_q_points), + rl2(n_q_points); + for (unsigned int k=0; k(&tmp_rhs[0], &rl1[0]); + eval_rhs.template values<1,false,false> (&rl1[0], &rl2[0]); + eval_rhs.template values<0,false,false> (&tmp_rhs[n_q_points], &rl1[0]); + eval_rhs.template gradients<1,false,true>(&rl1[0], &rl2[0]); + eval_rhs.template values<2,false,false> (&rl2[0], &tmp_rhs[0]); + eval_rhs.template values<0,false,false> (&tmp_rhs[2*n_q_points], &rl1[0]); + eval_rhs.template values<1,false,false> (&rl1[0], &rl2[0]); + eval_rhs.template gradients<2,false,true> (&rl2[0], &tmp_rhs[0]); + } + else if (dim == 2) + { + eval_rhs.template gradients<0,false,false>(&tmp_rhs[0], &rl1[0]); + eval_rhs.template values<1,false,false> (&rl1[0], &tmp_rhs[0]); + eval_rhs.template values<0,false,false> (&tmp_rhs[n_q_points], &rl1[0]); + eval_rhs.template gradients<1,false,true>(&rl1[0], &tmp_rhs[0]); + } + else + Assert(false, ExcNotImplemented()); + + tensor_product_matrix.apply_inverse(&rl1[0], &tmp_rhs[0]); + + for (unsigned int i=0; i