]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Finish the implementation of the criss-cross element.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 30 Jul 1998 09:06:37 +0000 (09:06 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 30 Jul 1998 09:06:37 +0000 (09:06 +0000)
git-svn-id: https://svn.dealii.org/trunk@462 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/source/fe/fe_lib.criss_cross.cc

index 1f6d31c68d73336b92ad552afcfcbfd66e47cfc3..b132df5a72e37435d620213769c7e5c2a5e321c9 100644 (file)
     od:
   od:
 
-  eq_sys := {(1-t)*x0 + t*x2 = (1-s)*x1 + s*x3,
-             (1-t)*y0 + t*y2 = (1-s)*y1 + s*y3}:
-  solution := solve (eq_sys, {s,t});
-
-  xs := subs (solution, (1-t)*x0 + t*x2):
-  ys := subs (solution, (1-t)*y0 + t*y2):
-  ps := array(1..2, [xs, ys]):
-
-  print ("writing data to files"):
-  readlib(C):
-  C(prolongation, filename=prolongation_2d):
-  C(ps, filename=crosspoint):
-
-  --------------------------------------------------------------------
-  
-  Postprocess the prolongation matrix by the commands
-  
-  perl -pi -e 's/\[(\d+)\]\[(\d+)\]\[(\d+)\]/[$1]($2,$3)/g;' prolongation_2d
-  perl -pi -e 's/.*= 0.0;\n//g;' prolongation_2d
-
------------------------------------------------------------------------------*/
-
-
-/*--------------------------------------
-    # these are the basis functions differentiated with respect to
+  # these are the basis functions differentiated with respect to
   # xi and eta. we need them for the computation of the jacobi
   # matrix, since we can't just differentiate a function.
   phi_xi[0] := proc(x,y) if(y<1-x) then -1;  else 0; fi; end:
   # define an array of the ansatz points in real space; the first
   # four are the vertices, the last one is the crossing point of
   # the two diagonals
-  x := array(0..4);
-  y := array(0..4);
+  print ("Computing cross point"):
+  x := array(0..4):
+  y := array(0..4):
 
   eq_sys := {(1-t)*x[0] + t*x[2] = (1-s)*x[1] + s*x[3],
              (1-t)*y[0] + t*y[2] = (1-s)*y[1] + s*y[3]}:
-  solution := solve (eq_sys, {s,t});
+  solution := solve (eq_sys, {s,t}):
 
   # set last point in dependence of the first four
   x[4] := subs (solution, (1-t)*x[0] + t*x[2]):
   # this is the mapping from the unit cell to the real cell, only for
   # completeness; we can't use it here, since phi[i] can't be
   # differentiated
-  x_real := sum(x[s]*phi[s], s=0..4):
-  y_real := sum(y[s]*phi[s], s=0..4):
+  x_real := simplify(sum(x[s]*phi[s], s=0..4)):
+  y_real := simplify(sum(y[s]*phi[s], s=0..4)):
 
   # correct form of the jacobi determinant:
   #   detJ :=   diff(x_real,xi)*diff(y_real,eta)
              detJ1(xi,eta) * detJ2(xi,eta) -
             detJ3(xi,eta) * detJ4(xi,eta);
           end:
-----------------------------------------------------------*/
+
+
+  # Now for the mass matrix: we divide the entire cell into four
+  # sectors:
+  #
+  # *-------*
+  # |\     /|
+  # | \ 3 / |
+  # |  \ /  |
+  # |4  *  2|
+  # |  / \  |
+  # | / 1 \ |
+  # |/     \|
+  # *-------*
+  #
+  # In each of these sectors, the Jacobi determinant is constant
+  # so that we can assemble the local mass matrix by summation
+  # over these four sectors. Since the basis functions are as of
+  # now only expressed as if-then-else statements, we have to
+  # express them for each sector separately and name them
+  # phi_s[i]. detJ_s denotes the Jacobi determinant on this sector.
+
+  print ("Computing mass matrix"):
+
+  mass_matrix := array (0..n_functions-1, 0..n_functions-1):
+  for i from 0 to n_functions-1 do
+    for j from 0 to n_functions-1 do
+      mass_matrix[i,j] := 0:
+    od:
+  od:
+  
+  # sector 1
+  phi_s[0] := 1-x-y:
+  phi_s[1] := x-y:
+  phi_s[2] := 0:
+  phi_s[3] := 0:
+  phi_s[4] := 1 - phi_s[0] - phi_s[1] - phi_s[2] - phi_s[3]:
+
+  detJ_s := simplify(detJ(1/2, 1/4)):
+  
+  for i from 0 to n_functions-1 do
+    for j from 0 to n_functions-1 do
+      # split integral over sector into the two parts
+      mass_matrix[i,j] := mass_matrix[i,j] +
+                          int(int(phi_s[i] * phi_s[j],
+                                  y=0..x),
+                             x=0..1/2) * detJ_s:
+      mass_matrix[i,j] := mass_matrix[i,j] +
+                          int(int(phi_s[i] * phi_s[j],
+                                  y=0..1-x),
+                             x=1/2..1) * detJ_s:
+    od:
+  od:  
+  
+  # sector 2
+  phi_s[0] := 0:
+  phi_s[1] := x-y:
+  phi_s[2] := x+y-1:
+  phi_s[3] := 0:
+  phi_s[4] := 1 - phi_s[0] - phi_s[1] - phi_s[2] - phi_s[3]:
+
+  detJ_s := simplify(detJ(3/4, 1/2)):
+  
+  for i from 0 to n_functions-1 do
+    for j from 0 to n_functions-1 do
+      # split integral over sector into the two parts
+      mass_matrix[i,j] := mass_matrix[i,j] +
+                          int(int(phi_s[i] * phi_s[j],
+                                  y=0..x),
+                             x=0..1/2) * detJ_s:
+      mass_matrix[i,j] := mass_matrix[i,j] +
+                          int(int(phi_s[i] * phi_s[j],
+                                  y=0..1-x),
+                             x=1/2..1) * detJ_s:
+    od:
+  od:  
+
+  # sector 3
+  phi_s[0] := 0:
+  phi_s[1] := 0:
+  phi_s[2] := x+y-1:
+  phi_s[3] := y-x:
+  phi_s[4] := 1 - phi_s[0] - phi_s[1] - phi_s[2] - phi_s[3]:
+
+  detJ_s := simplify(detJ(1/2, 3/4)):
+  
+  for i from 0 to n_functions-1 do
+    for j from 0 to n_functions-1 do
+      # split integral over sector into the two parts
+      mass_matrix[i,j] := mass_matrix[i,j] +
+                          int(int(phi_s[i] * phi_s[j],
+                                  y=0..x),
+                             x=0..1/2) * detJ_s:
+      mass_matrix[i,j] := mass_matrix[i,j] +
+                          int(int(phi_s[i] * phi_s[j],
+                                  y=0..1-x),
+                             x=1/2..1) * detJ_s:
+    od:
+  od:  
+
+  # sector 4
+  phi_s[0] := 1-x-y:
+  phi_s[1] := 0:
+  phi_s[2] := 0:
+  phi_s[3] := y-x:
+  phi_s[4] := 1 - phi_s[0] - phi_s[1] - phi_s[2] - phi_s[3]:
+
+  detJ_s := simplify(detJ(1/4, 1/2)):
+  
+  for i from 0 to n_functions-1 do
+    for j from 0 to n_functions-1 do
+      # split integral over sector into the two parts
+      mass_matrix[i,j] := mass_matrix[i,j] +
+                          int(int(phi_s[i] * phi_s[j],
+                                  y=0..x),
+                             x=0..1/2) * detJ_s:
+      mass_matrix[i,j] := mass_matrix[i,j] +
+                          int(int(phi_s[i] * phi_s[j],
+                                  y=0..1-x),
+                             x=1/2..1) * detJ_s:
+    od:
+  od:  
+  
+  print ("writing data to files"):
+  readlib(C):
+  C(prolongation, filename=prolongation_2d):
+  C(array(1..2, [x[4], y[4]]), optimized, filename=crosspoint_2d):
+  C(mass_matrix, optimized, filename=massmatrix_2d):
+  
+  --------------------------------------------------------------------
+  
+  Postprocess the files by the commands
+  
+  perl -pi -e 's/\[(\d+)\]\[(\d+)\]\[(\d+)\]/[$1]($2,$3)/g;' prolongation_2d
+  perl -pi -e 's/.*= 0.0;\n//g;' prolongation_2d
+  perl -pi -e 's/\[(\d+)\]\[(\d+)\]/($1,$2)/g;' massmatrix_2d
+  perl -pi -e 's/(t\d+) =/const double $1 =/g;' massmatrix_2d
+
+-----------------------------------------------------------------------------*/
+
+
+
 
 
 
@@ -451,22 +569,144 @@ void FECrissCross<2>::get_face_ansatz_points (const DoFHandler<2>::face_iterator
 
 
 template <>
-void FECrissCross<2>::get_local_mass_matrix (const DoFHandler<2>::cell_iterator &,
+void FECrissCross<2>::get_local_mass_matrix (const DoFHandler<2>::cell_iterator &cell,
                                             const Boundary<2> &,
-                                            dFMatrix &local_mass_matrix) const {
-  Assert (local_mass_matrix.n() == total_dofs,
-         ExcWrongFieldDimension(local_mass_matrix.n(),total_dofs));
-  Assert (local_mass_matrix.m() == total_dofs,
-         ExcWrongFieldDimension(local_mass_matrix.m(),total_dofs));
-
-                                  // in this special element, some of the
-                                  // entries are zero (which is not the
-                                  // case for most other elements, so
-                                  // we first reset all elements and only
-                                  // fill in those that are nonzero
-  local_mass_matrix.clear ();
-  
-  Assert (false, ExcNotUseful());
+                                            dFMatrix &mass_matrix) const {
+  Assert (mass_matrix.n() == total_dofs,
+         ExcWrongFieldDimension(mass_matrix.n(),total_dofs));
+  Assert (mass_matrix.m() == total_dofs,
+         ExcWrongFieldDimension(mass_matrix.m(),total_dofs));
+
+  const double x[4] = { cell->vertex(0)(0),
+                       cell->vertex(1)(0),
+                       cell->vertex(2)(0),
+                       cell->vertex(3)(0)  };
+  const double y[4] = { cell->vertex(0)(1),
+                       cell->vertex(1)(1),
+                       cell->vertex(2)(1),
+                       cell->vertex(3)(1)  };
+
+  const double t1 = x[3]*x[2];
+  const double t2 = y[1]*y[1];
+  const double t5 = x[0]*x[0];
+  const double t7 = t5*y[3]*y[2];
+  const double t8 = x[3]*x[1];
+  const double t9 = y[0]*y[0];
+  const double t10 = t8*t9;
+  const double t11 = t5*y[1];
+  const double t12 = t11*y[3];
+  const double t13 = t11*y[2];
+  const double t14 = x[1]*x[1];
+  const double t15 = t14*y[0];
+  const double t19 = x[2]*x[1]*t9;
+  const double t20 = t1*t9;
+  const double t24 = x[0]*x[2];
+  const double t25 = t24*t2;
+  const double t29 = t15*y[2];
+  const double t30 = x[3]*y[0];
+  const double t31 = x[1]*y[2];
+  const double t32 = t30*t31;
+  const double t33 = x[2]*y[3];
+  const double t34 = x[1]*y[0];
+  const double t35 = t33*t34;
+  const double t37 = x[2]*y[1];
+  const double t38 = t37*t34;
+  const double t39 = y[1]*x[1];
+  const double t42 = -2.0*t1*t2-t7+t10+t12+t13+2.0*t15*y[3]+t19-t20-2.0*t14*y[2]*y[3]+
+                    t25+2.0*x[0]*t2*x[3]+t29-t32-2.0*t35-t38-2.0*t39*t30;
+  const double t43 = y[3]*x[1];
+  const double t46 = x[3]*y[2];
+  const double t49 = t37*t30;
+  const double t51 = x[0]*y[1];
+  const double t52 = t51*t46;
+  const double t54 = t51*t33;
+  const double t55 = x[0]*y[0];
+  const double t56 = t55*t37;
+  const double t59 = t51*t30;
+  const double t60 = t51*t31;
+  const double t61 = t55*t31;
+  const double t62 = t55*t46;
+  const double t63 = x[0]*y[3];
+  const double t64 = t63*t34;
+  const double t65 = t55*t33;
+  const double t66 = t63*t31;
+  const double t72 = 2.0*t37*t43+2.0*t39*t46+3.0*t49-2.0*t52-t54-t56-2.0*t51*t43-t59-t60
+                    -t61+t62-t64+t65+3.0*t66-t14*t9-t5*t2+2.0*t51*t34;
+  const double t75 = 1/(t51-t63-t37+t33-t34+t31+t30-t46);
+  const double t76 = (t42+t72)*t75;
+  const double t77 = y[3]*y[3];
+  const double t81 = x[3]*x[3];
+  const double t82 = y[1]*t81;
+  const double t86 = t81*y[0]*y[2];
+  const double t90 = t24*t77;
+  const double t94 = -t7-t10-t12-2.0*x[0]*t77*x[1]+t13+t19+2.0*t82*y[2]-t86-t20+t81*t9
+                    -2.0*t82*y[0]-t90+t32-3.0*t35+2.0*t49-3.0*t52;
+  const double t96 = t63*t46;
+  const double t97 = t30*t33;
+  const double t98 = x[3]*y[3];
+  const double t114 = t54-t56+t59-t61+t62+t64+t65+2.0*t66+t96+t97+2.0*t51*t98-2.0*t63*
+                     t30-2.0*y[1]*x[3]*t33+2.0*t98*t34-2.0*t98*t31+2.0*t77*x[1]*x[2]+t5*t77;
+  const double t116 = (t94+t114)*t75;
+  const double t118 = t76/24.0;
+  const double t119 = t116/24.0;
+  const double t121 = -t118+t116/8.0;
+  const double t122 = x[0]*y[2];
+  const double t123 = t122*t37;
+  const double t124 = t122*t33;
+  const double t125 = y[2]*y[2];
+  const double t126 = x[0]*t125;
+  const double t127 = t126*x[1];
+  const double t128 = t126*x[3];
+  const double t129 = x[2]*x[2];
+  const double t131 = y[1]*t129;
+  const double t132 = t131*y[0];
+  const double t133 = t131*y[3];
+  const double t134 = -t25+t60+t123+t54-t52-t124-t127+t128+t129*t2-t132-t133;
+  const double t135 = t37*t46;
+  const double t139 = x[2]*y[0];
+  const double t140 = t139*t46;
+  const double t141 = t8*t125;
+  const double t143 = t129*y[0]*y[3];
+  const double t144 = t33*t31;
+  const double t145 = t139*t31;
+  const double t146 = t38+t135-2.0*t37*t31-t29+t14*t125-t140-t35-t141+t143+t144+t32+t145
+                     ;
+  const double t148 = (t134+t146)*t75;
+  const double t150 = t148/24.0;
+  const double t152 = -t118+t148/8.0;
+  const double t153 = t123-t54+t66+t128-t96-t124-t127+t90+t49-t132-t135;
+  const double t158 = t133+2.0*t33*t46+t86-t81*t125-t129*t77+t141-t97-t32-t144-t140+t143
+                     +t145;
+  const double t160 = (t153+t158)*t75;
+  const double t162 = t160/24.0;
+  const double t164 = 7.0/24.0*t160;
+  const double t165 = -5.0/24.0*t148+t164;
+  const double t168 = t164-5.0/24.0*t116;
+  mass_matrix(0,0) = -t76/12.0+t116/12.0;
+  mass_matrix(0,1) = -t118;
+  mass_matrix(0,2) = 0.0;
+  mass_matrix(0,3) = -t119;
+  mass_matrix(0,4) = t121;
+  mass_matrix(1,0) = -t118;
+  mass_matrix(1,1) = -t76/12.0+t148/12.0;
+  mass_matrix(1,2) = -t150;
+  mass_matrix(1,3) = 0.0;
+  mass_matrix(1,4) = t152;
+  mass_matrix(2,0) = 0.0;
+  mass_matrix(2,1) = -t150;
+  mass_matrix(2,2) = t148/12.0-t160/12.0;
+  mass_matrix(2,3) = -t162;
+  mass_matrix(2,4) = t165;
+  mass_matrix(3,0) = -t119;
+  mass_matrix(3,1) = 0.0;
+  mass_matrix(3,2) = -t162;
+  mass_matrix(3,3) = -t160/12.0+t116/12.0;
+  mass_matrix(3,4) = t168;
+  mass_matrix(4,0) = t121;
+  mass_matrix(4,1) = t152;
+  mass_matrix(4,2) = t165;
+  mass_matrix(4,3) = t168;
+  mass_matrix(4,4) = -t76/12.0+7.0/12.0*t148-17.0/12.0*t160+7.0/12.0*t116;
 };
 
 

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.