/*------------------------------- 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,
template <>
void FiniteElement<1>::get_support_points (const DoFHandler<1>::cell_iterator &cell,
- const Boundary<1> &,
- vector<Point<1> > &support_points) const {
+ const Boundary<1> &,
+ vector<Point<1> > &support_points) const {
Assert (support_points.size() == total_dofs,
ExcWrongFieldDimension(support_points.size(), total_dofs));
// compute support points. The first ones
#include <algorithm>
-/*--------------------------------- 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
-*/
-
-
+template <>
+Tensor<2,1>
+FECubicSub<1>::shape_grad_grad (const unsigned int i,
+ const Point<1> &p) const
+{
+ Assert (i<total_dofs, ExcInvalidIndex(i));
+
+ const double xi = p(0);
+ Tensor<2,1> 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<Point<1> > &unit_points) const {
FiniteElement<1>::get_unit_support_points (unit_points);
+template <>
+Tensor<2,1>
+FEQuadraticSub<1>::shape_grad_grad (const unsigned int i,
+ const Point<1> &) const
+{
+ Assert((i<total_dofs), ExcInvalidIndex(i));
+
+ Tensor<2,1> 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<Point<1> > &unit_points) const {
FiniteElement<1>::get_unit_support_points (unit_points);
#include <algorithm>
-/*--------------------------------- 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
-*/
-
-
-
-
+template <>
+Tensor<2,1>
+FEQuarticSub<1>::shape_grad_grad (const unsigned int i,
+ const Point<1> &p) const
+{
+ Assert (i<total_dofs, ExcInvalidIndex(i));
+
+ const double xi = p(0);
+ Tensor<2,1> 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<Point<1> > &unit_points) const {
FiniteElement<1>::get_unit_support_points (unit_points);