From: wolf Date: Fri, 6 Nov 1998 16:16:32 +0000 (+0000) Subject: Implement second derivatives for 1D. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a7d44d800fa4f672de96c5dd90868a7700c4b431;p=dealii-svn.git Implement second derivatives for 1D. git-svn-id: https://svn.dealii.org/trunk@649 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index 1726959110..f3379ebea6 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -162,11 +162,12 @@ bool FiniteElementBase::operator == (const FiniteElementBase &f) const /*------------------------------- FiniteElement ----------------------*/ +#if deal_II_dimension == 1 + // declare this function to be explicitely specialized before first use // egcs wants this, but gcc2.8.1 produces an internal compiler error, so // we drop this declaration again for the time being -#if deal_II_dimension == 1 //template <> //void FiniteElement<1>::get_support_points (const DoFHandler<1>::cell_iterator &cell, @@ -286,8 +287,8 @@ void FiniteElement<1>::get_unit_support_points (vector > &support_point template <> void FiniteElement<1>::get_support_points (const DoFHandler<1>::cell_iterator &cell, - const Boundary<1> &, - vector > &support_points) const { + const Boundary<1> &, + vector > &support_points) const { Assert (support_points.size() == total_dofs, ExcWrongFieldDimension(support_points.size(), total_dofs)); // compute support points. The first ones diff --git a/deal.II/deal.II/source/fe/fe_lib.cubic.cc b/deal.II/deal.II/source/fe/fe_lib.cubic.cc index 654489f5f0..1216637d63 100644 --- a/deal.II/deal.II/source/fe/fe_lib.cubic.cc +++ b/deal.II/deal.II/source/fe/fe_lib.cubic.cc @@ -9,120 +9,6 @@ #include -/*--------------------------------- For 1d --------------------------------- - -- Use the following maple script to generate the basis functions, - -- gradients and prolongation matrices as well as the mass matrix. - -- Make sure that the files do not exists beforehand, since output - -- is appended instead of overwriting previous contents. - -- - -- You should only have to change the very first lines for polynomials - -- of higher order. - -------------------------------------------------------------------------- - n_functions := 4; - - support_points := array(0..n_functions-1); - support_points[0] := 0; - support_points[1] := 1; - support_points[2] := 1/3; - support_points[3] := 2/3; - - phi_polynom := array(0..n_functions-1); - grad_phi_polynom := array(0..n_functions-1); - local_mass_matrix := array(0..n_functions-1, 0..n_functions-1); - - for i from 0 to n_functions-1 do - # note that the interp function wants vectors indexed from - # one and not from zero. - values := array(1..n_functions); - for j from 1 to n_functions do - values[j] := 0; - od; - values[i+1] := 1; - - shifted_support_points := array (1..n_functions); - for j from 1 to n_functions do - shifted_support_points[j] := support_points[j-1]; - od; - - phi_polynom[i] := interp (shifted_support_points, values, xi); - grad_phi_polynom[i] := diff(phi_polynom[i], xi); - od; - - phi:= proc(i,x) subs(xi=x, phi_polynom[i]); end; - - - points[0] := array(0..n_functions-1); - points[1] := array(0..n_functions-1); - for i from 0 to n_functions-1 do - points[0][i] := support_points[i]/2; - points[1][i] := support_points[i]/2+1/2; - od; - - prolongation := array(0..1,0..n_functions-1, 0..n_functions-1); - - for i from 0 to 1 do - for j from 0 to n_functions-1 do - for k from 0 to n_functions-1 do - prolongation[i,j,k] := phi(k, points[i][j]); - od; - od; - od; - - - # to get the restriction (interpolation) matrices, evaluate - # the basis functions on the child cells at the global - # interpolation points - child_phi[0] := proc(i, point) - if ((point<0) or (point>1/2)) then - 0: - else - phi(i,2*point): - fi: - end: - child_phi[1] := proc(i, point) - if ((point<1/2) or (point>1)) then - 0: - else - phi(i,2*point-1): - fi: - end: - restriction := array(0..1,0..n_functions-1, 0..n_functions-1); - for child from 0 to 1 do - for j from 0 to n_functions-1 do - for k from 0 to n_functions-1 do - restriction[child,j,k] := child_phi[child](k, support_points[j]): - od: - od: - od: - - - for i from 0 to n_functions-1 do - for j from 0 to n_functions-1 do - local_mass_matrix[i,j] := int(phi_polynom[i] * phi_polynom[j] * h, - xi=0..1); - od; - od; - - readlib(C); - C(phi_polynom, filename=shape_value_1d); - C(grad_phi_polynom, filename=shape_grad_1d); - C(prolongation, filename=prolongation_1d); - C(restriction, filename=restriction_1d); - C(local_mass_matrix, optimized, filename=massmatrix_1d); - - ----------------------------------------------------------------------- - Use the following perl scripts to convert the output into a - suitable format: - - perl -pi -e 's/phi_polynom\[(\d)\] =/case $1: return/g;' shape_value_1d - perl -pi -e 's/grad_phi_polynom\[(\d)\] = (.*);/case $1: return Point<1>($2);/g;' shape_grad_1d - perl -pi -e 's/\[(\d+)\]\[(\d)\]/($1,$2)/g;' massmatrix_1d - perl -pi -e 's/\[(\d+)\]\[(\d)\]\[(\d)\]/[$1]($2,$3)/g;' prolongation_1d - perl -pi -e 's/\[(\d+)\]\[(\d)\]\[(\d)\]/[$1]($2,$3)/g;' restriction_1d - perl -pi -e 's/(t\d+) =/const double $1 =/g;' massmatrix_1d -*/ - - @@ -239,6 +125,28 @@ FECubicSub<1>::shape_grad (const unsigned int i, +template <> +Tensor<2,1> +FECubicSub<1>::shape_grad_grad (const unsigned int i, + const Point<1> &p) const +{ + Assert (i return_value; + switch (i) + { + case 0: return_value[0][0] = -27.0*xi+18.0; + case 1: return_value[0][0] = 27.0*xi-9.0; + case 2: return_value[0][0] = 81.0*xi-45.0; + case 3: return_value[0][0] = -81.0*xi+36.0; + }; + + return return_value; +}; + + + template <> void FECubicSub<1>::get_unit_support_points (vector > &unit_points) const { FiniteElement<1>::get_unit_support_points (unit_points); diff --git a/deal.II/deal.II/source/fe/fe_lib.quadratic.cc b/deal.II/deal.II/source/fe/fe_lib.quadratic.cc index 5ace45a8ef..49236de6d6 100644 --- a/deal.II/deal.II/source/fe/fe_lib.quadratic.cc +++ b/deal.II/deal.II/source/fe/fe_lib.quadratic.cc @@ -135,6 +135,31 @@ FEQuadraticSub<1>::shape_grad(const unsigned int i, +template <> +Tensor<2,1> +FEQuadraticSub<1>::shape_grad_grad (const unsigned int i, + const Point<1> &) const +{ + Assert((i return_value; + switch (i) + { + case 0: + return_value[0][0] = 4; + break; + case 1: + return_value[0][0] = 4; + break; + case 2: + return_value[0][0] = -8; + break; + } + return return_value; +}; + + + template <> void FEQuadraticSub<1>::get_unit_support_points (vector > &unit_points) const { FiniteElement<1>::get_unit_support_points (unit_points); diff --git a/deal.II/deal.II/source/fe/fe_lib.quartic.cc b/deal.II/deal.II/source/fe/fe_lib.quartic.cc index 134f9d9faa..7cd303f6e1 100644 --- a/deal.II/deal.II/source/fe/fe_lib.quartic.cc +++ b/deal.II/deal.II/source/fe/fe_lib.quartic.cc @@ -9,126 +9,6 @@ #include -/*--------------------------------- For 1d --------------------------------- - -- Use the following maple script to generate the basis functions, - -- gradients and prolongation matrices as well as the mass matrix. - -- Make sure that the files do not exists beforehand, since output - -- is appended instead of overwriting previous contents. - -- - -- You should only have to change the very first lines for polynomials - -- of higher order. - -------------------------------------------------------------------------- - n_functions := 5; - - support_points := array(0..n_functions-1); - support_points[0] := 0; - support_points[1] := 1; - support_points[2] := 1/4; - support_points[3] := 2/4; - support_points[4] := 3/4; - - phi_polynom := array(0..n_functions-1); - grad_phi_polynom := array(0..n_functions-1); - local_mass_matrix := array(0..n_functions-1, 0..n_functions-1); - - for i from 0 to n_functions-1 do - # note that the interp function wants vectors indexed from - # one and not from zero. - values := array(1..n_functions); - for j from 1 to n_functions do - values[j] := 0; - od; - values[i+1] := 1; - - shifted_support_points := array (1..n_functions); - for j from 1 to n_functions do - shifted_support_points[j] := support_points[j-1]; - od; - - phi_polynom[i] := interp (shifted_support_points, values, xi); - grad_phi_polynom[i] := diff(phi_polynom[i], xi); - od; - - phi:= proc(i,x) subs(xi=x, phi_polynom[i]); end; - - - points[0] := array(0..n_functions-1); - points[1] := array(0..n_functions-1); - for i from 0 to n_functions-1 do - points[0][i] := support_points[i]/2; - points[1][i] := support_points[i]/2+1/2; - od; - - prolongation := array(0..1,0..n_functions-1, 0..n_functions-1); - - for i from 0 to 1 do - for j from 0 to n_functions-1 do - for k from 0 to n_functions-1 do - prolongation[i,j,k] := phi(k, points[i][j]); - od; - od; - od; - - for i from 0 to n_functions-1 do - for j from 0 to n_functions-1 do - local_mass_matrix[i,j] := int(phi_polynom[i] * phi_polynom[j] * h, - xi=0..1); - od; - od; - - # to get the restriction (interpolation) matrices, evaluate - # the basis functions on the child cells at the global - # interpolation points - child_phi[0] := proc(i, point) - if ((point<0) or (point>1/2)) then - 0: - else - phi(i,2*point): - fi: - end: - child_phi[1] := proc(i, point) - if ((point<1/2) or (point>1)) then - 0: - else - phi(i,2*point-1): - fi: - end: - restriction := array(0..1,0..n_functions-1, 0..n_functions-1); - for child from 0 to 1 do - for j from 0 to n_functions-1 do - for k from 0 to n_functions-1 do - restriction[child,j,k] := child_phi[child](k, support_points[j]): - od: - od: - od: - - - readlib(C); - C(phi_polynom, filename=shape_value_1d); - C(grad_phi_polynom, filename=shape_grad_1d); - C(prolongation, filename=prolongation_1d); - C(restriction, filename=restriction_1d); - C(local_mass_matrix, optimized, filename=massmatrix_1d); - - ----------------------------------------------------------------------- - Use the following perl scripts to convert the output into a - suitable format: - - perl -pi -e 's/([^;])\n/$1/g;' shape_value_1d - perl -pi -e 's/([^;])\n/$1/g;' shape_grad_1d - perl -pi -e 's/phi_polynom\[(\d)\] =/case $1: return/g;' shape_value_1d - perl -pi -e 's/grad_phi_polynom\[(\d)\] = (.*);/case $1: return Point<1>($2);/g;' shape_grad_1d - perl -pi -e 's/\[(\d+)\]\[(\d)\]/($1,$2)/g;' massmatrix_1d - perl -pi -e 's/\[(\d+)\]\[(\d)\]\[(\d)\]/[$1]($2,$3)/g;' prolongation_1d - perl -pi -e 's/.*= 0.0;\n//g;' prolongation_1d - perl -pi -e 's/\[(\d+)\]\[(\d)\]\[(\d)\]/[$1]($2,$3)/g;' restriction_1d - perl -pi -e 's/.*= 0.0;\n//g;' restriction_1d - perl -pi -e 's/(t\d+) =/const double $1 =/g;' massmatrix_1d -*/ - - - - @@ -214,6 +94,29 @@ FEQuarticSub<1>::shape_grad(const unsigned int i, +template <> +Tensor<2,1> +FEQuarticSub<1>::shape_grad_grad (const unsigned int i, + const Point<1> &p) const +{ + Assert (i return_value; + switch (i) + { + case 0: return_value[0][0] = 128.0*xi*xi-160.0*xi+140.0/3.0; + case 1: return_value[0][0] = 128.0*xi*xi-96.0*xi+44.0/3.0; + case 2: return_value[0][0] = -512.0*xi*xi+576.0*xi-416.0/3.0; + case 3: return_value[0][0] = 768.0*xi*xi-768.0*xi+152.0; + case 4: return_value[0][0] = -512.0*xi*xi+448.0*xi-224.0/3.0; + }; + + return return_value; +}; + + + template <> void FEQuarticSub<1>::get_unit_support_points (vector > &unit_points) const { FiniteElement<1>::get_unit_support_points (unit_points);