]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Implement second derivatives for 1D.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 6 Nov 1998 16:16:32 +0000 (16:16 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 6 Nov 1998 16:16:32 +0000 (16:16 +0000)
git-svn-id: https://svn.dealii.org/trunk@649 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/source/fe/fe.cc
deal.II/deal.II/source/fe/fe_lib.cubic.cc
deal.II/deal.II/source/fe/fe_lib.quadratic.cc
deal.II/deal.II/source/fe/fe_lib.quartic.cc

index 1726959110ae816614ab552d2ab4aecf9b5bdf13..f3379ebea65f98bde50e9f1fab48c487d1b80c92 100644 (file)
@@ -162,11 +162,12 @@ bool FiniteElementBase<dim>::operator == (const FiniteElementBase<dim> &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<Point<1> > &support_point
 
 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
index 654489f5f05455c4364372f52b35127349e5243b..1216637d6343a7c6f9c0ad07e1c86c54b4a868d8 100644 (file)
@@ -9,120 +9,6 @@
 #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
-*/
-
-
 
 
 
@@ -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<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);
index 5ace45a8efd83f27c4d0e6c34892369ade0ccc02..49236de6d6c4e71e0aa197055ae2ffb3edbc20ab 100644 (file)
@@ -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<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);
index 134f9d9faa321c1a595f7910d63540ced4de6a31..7cd303f6e1bd613dd241fb8ea48239b80f963162 100644 (file)
@@ -9,126 +9,6 @@
 #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
-*/
-
-
-
-
 
 
 
@@ -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<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);

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.