# Copyright W. Bangerth, University of Heidelberg, 1998
# Global options for all compilations (compilations in /source now
-# use other options!)
+# use other options! THEY REALLY SHOULD NOT (gk))
+
+# This file uses the variable $D, which is the deal home directory,
+# that is, the parent directory of this directory.
+#
+# Please set $D to the right path before including this file!
deal_II_dimension=2
+vpath %.a: $D/base/lib
+vpath %.a: $D/lac/lib
+vpath %.a: $D/deal.II/lib
+
CXX = c++
-INCLUDE = -I../../include -I../../../lac/include \
- -I../../../base/include -I../../../mia/include
+INCLUDE = -I$D/deal.II/include -I$D/lac/include \
+ -I$D/base/include
CXXFLAGS.g= -DDEBUG -g -Wall -W -pedantic -Wconversion \
-Winline -Woverloaded-virtual -fno-rtti \
$(INCLUDE) -Ddeal_II_dimension=$(deal_II_dimension)
CXX = /home/wolf/bin/gcc/bin/c++
endif
-ifeq ($(shell uname),SunOS)
-CXX = /usr/local/gcc/egcs-19980803/bin/c++
-endif
-
-
%.go : %.cc Makefile
- @echo ============================ Compiling with debugging information: $<
+ @echo ======================== Compiling (DEBUG): $<
# @echo $(CXX) ... -c $< -o $@
@$(CXX) $(CXXFLAGS.g) -c $< -o $@
%.o : %.cc Makefile
- @echo ============================ Compiling with optimization: $<
+ @echo ======================== Compiling (OPT): $<
# @echo $(CXX) ... -c $< -o $@
@$(CXX) $(CXXFLAGS) -c $< -o $@
source-clean:
cd source ; make clean
-
+TAGS:
+ etags include/*/* source/*/*
.PHONY: all deal.II doc examples
* The question what a degree of freedom on the boundary is, is not so easy.
* It should really be a degree of freedom of which the respective basis
* function has nonzero values on the boundary. At least for Lagrange elements
- * this definition is equal to the statement that the off-point of the ansatz
+ * this definition is equal to the statement that the off-point of the trial
* function, i.e. the point where the function assumes its nominal value (for
* Lagrange elements this is the point where it has the function value #1#), is
* located on the boundary. We do not check this directly, the criterion is
* not on the boundary, the value of
* #mapping[dof]# is #-1#. This function is
* mainly used when setting up matrices and
- * vectors on the boundary from the ansatz
+ * vectors on the boundary from the trial
* functions, which have global numbers,
* while the matrices and vectors use
- * numbers of the ansatz functions local
+ * numbers of the trial functions local
* to the boundary.
*
* Prior content of #mapping# is deleted.
* linear on each triangles and vanishes at the faces of the quadrilateral.
*
* Now, on the unit element, these basis functions are the same as for
- * a triangular ansatz space, namely the class of ${\cal P}_1$ Lagrange
+ * a triangular trial space, namely the class of ${\cal P}_1$ Lagrange
* elements. Due to the arrangement of the four triangles on the
* quadrilateral, it is clear that considering the whole triangulation
* of the domain, always four triangles meet at the points which
* between unit and real cell, the mapping of each of the four triangles
* is linear (note that this does not hold for the whole element; the
* mapping of the quadrilaterals is only piecewise linear and globally
- * continuous). It also proves that the ansatz space using this element
- * is equivalent to the ansatz space using triangles and linear elements.
+ * continuous). It also proves that the trial space using this element
+ * is equivalent to the trial space using triangles and linear elements.
*
* Since in one space dimension, this element equals two linear elements,
* i.e. a linear ansatz on a mesh once more refined than the present one,
* Refer to the base class for detailed
* information on this function.
*/
- virtual void get_unit_ansatz_points (vector<Point<dim> > &ansatz_points) const;
+ virtual void get_unit_support_points (vector<Point<dim> > &support_points) const;
/**
* Refer to the base class for detailed
* information on this function.
*/
- virtual void get_ansatz_points (const DoFHandler<dim>::cell_iterator &cell,
+ virtual void get_support_points (const DoFHandler<dim>::cell_iterator &cell,
const Boundary<dim> &boundary,
- vector<Point<dim> > &ansatz_points) const;
+ vector<Point<dim> > &support_points) const;
/**
* Refer to the base class for detailed
* information on this function.
*/
- virtual void get_face_ansatz_points (const DoFHandler<dim>::face_iterator &face,
+ virtual void get_face_support_points (const DoFHandler<dim>::face_iterator &face,
const Boundary<dim> &boundary,
- vector<Point<dim> > &ansatz_points) const;
+ vector<Point<dim> > &support_points) const;
/**
* Refer to the base class for detailed
const vector<Point<dim> > &unit_points,
vector<dFMatrix> &jacobians,
const bool compute_jacobians,
- vector<Point<dim> > &ansatz_points,
- const bool compute_ansatz_points,
+ vector<Point<dim> > &support_points,
+ const bool compute_support_points,
vector<Point<dim> > &q_points,
const bool compute_q_points,
const dFMatrix &shape_values_transform,
* Refer to the base class for detailed
* information on this function.
*/
- virtual void get_unit_ansatz_points (vector<Point<dim> > &ansatz_points) const;
+ virtual void get_unit_support_points (vector<Point<dim> > &support_points) const;
/**
* Refer to the base class for detailed
* information on this function.
*/
- virtual void get_ansatz_points (const DoFHandler<dim>::cell_iterator &cell,
+ virtual void get_support_points (const DoFHandler<dim>::cell_iterator &cell,
const Boundary<dim> &boundary,
- vector<Point<dim> > &ansatz_points) const;
+ vector<Point<dim> > &support_points) const;
/**
* Refer to the base class for detailed
* information on this function.
*/
- virtual void get_face_ansatz_points (const DoFHandler<dim>::face_iterator &face,
+ virtual void get_face_support_points (const DoFHandler<dim>::face_iterator &face,
const Boundary<dim> &boundary,
- vector<Point<dim> > &ansatz_points) const;
+ vector<Point<dim> > &support_points) const;
/**
* Refer to the base class for detailed
* Refer to the base class for detailed
* information on this function.
*/
- virtual void get_unit_ansatz_points (vector<Point<dim> > &ansatz_points) const;
+ virtual void get_unit_support_points (vector<Point<dim> > &support_points) const;
/**
* Refer to the base class for detailed
* information on this function.
*/
- virtual void get_ansatz_points (const DoFHandler<dim>::cell_iterator &cell,
+ virtual void get_support_points (const DoFHandler<dim>::cell_iterator &cell,
const Boundary<dim> &boundary,
- vector<Point<dim> > &ansatz_points) const;
+ vector<Point<dim> > &support_points) const;
/**
* Refer to the base class for detailed
* information on this function.
*/
- virtual void get_face_ansatz_points (const DoFHandler<dim>::face_iterator &face,
+ virtual void get_face_support_points (const DoFHandler<dim>::face_iterator &face,
const Boundary<dim> &boundary,
- vector<Point<dim> > &ansatz_points) const;
+ vector<Point<dim> > &support_points) const;
/**
* Refer to the base class for detailed
* Refer to the base class for detailed
* information on this function.
*/
- virtual void get_unit_ansatz_points (vector<Point<dim> > &ansatz_points) const;
+ virtual void get_unit_support_points (vector<Point<dim> > &support_points) const;
/**
* Refer to the base class for detailed
* information on this function.
*/
- virtual void get_ansatz_points (const DoFHandler<dim>::cell_iterator &cell,
+ virtual void get_support_points (const DoFHandler<dim>::cell_iterator &cell,
const Boundary<dim> &boundary,
- vector<Point<dim> > &ansatz_points) const;
+ vector<Point<dim> > &support_points) const;
/**
* Refer to the base class for detailed
* information on this function.
*/
- virtual void get_face_ansatz_points (const DoFHandler<dim>::face_iterator &face,
+ virtual void get_face_support_points (const DoFHandler<dim>::face_iterator &face,
const Boundary<dim> &boundary,
- vector<Point<dim> > &ansatz_points) const;
+ vector<Point<dim> > &support_points) const;
/**
* Refer to the base class for detailed
* Refer to the base class for detailed
* information on this function.
*/
- virtual void get_unit_ansatz_points (vector<Point<dim> > &ansatz_points) const;
+ virtual void get_unit_support_points (vector<Point<dim> > &support_points) const;
/**
* Refer to the base class for detailed
* information on this function.
*/
- virtual void get_ansatz_points (const DoFHandler<dim>::cell_iterator &cell,
+ virtual void get_support_points (const DoFHandler<dim>::cell_iterator &cell,
const Boundary<dim> &boundary,
- vector<Point<dim> > &ansatz_points) const;
+ vector<Point<dim> > &support_points) const;
/**
* Refer to the base class for detailed
* information on this function.
*/
- virtual void get_face_ansatz_points (const DoFHandler<dim>::face_iterator &face,
+ virtual void get_face_support_points (const DoFHandler<dim>::face_iterator &face,
const Boundary<dim> &boundary,
- vector<Point<dim> > &ansatz_points) const;
+ vector<Point<dim> > &support_points) const;
/**
* Refer to the base class for detailed
const vector<Point<dim> > &unit_points,
vector<dFMatrix> &jacobians,
const bool compute_jacobians,
- vector<Point<dim> > &ansatz_points,
- const bool compute_ansatz_points,
+ vector<Point<dim> > &support_points,
+ const bool compute_support_points,
vector<Point<dim> > &q_points,
const bool compute_q_points,
const dFMatrix &shape_values_transform,
* It can be understood as a gradient recovery estimator; see the survey
* of Ainsworth for a complete discussion.
*
- * It seem as if this error estimator should only be valid for linear ansatz
- * spaces, and there are indications that for higher order ansatz spaces the
+ * It seem as if this error estimator should only be valid for linear trial
+ * spaces, and there are indications that for higher order trial spaces the
* integrals computed here show superconvergence properties, i.e. they tend
* to zero faster than the error itself, thus ruling out the values as error
* indicators.
* square of the jumps) of each cell and take the square root.
*
* The integration is done using a quadrature formula on the face.
- * For linear ansatz functions (#FELinear#), the #Gauss2# or even the
+ * For linear trial functions (#FELinear#), the #Gauss2# or even the
* #Midpoint# rule will suffice. For higher order elements, it is
* necessary to utilize higher order quadrature formulae as well.
*
n_functions := 5:
- # note: ansatz_points[i] is a vector which is indexed from
+ # note: support_points[i] is a vector which is indexed from
# one and not from zero!
- ansatz_points[0] := [0,0]:
- ansatz_points[1] := [1,0]:
- ansatz_points[2] := [1,1]:
- ansatz_points[3] := [0,1]:
- ansatz_points[4] := [1/2,1/2]:
+ support_points[0] := [0,0]:
+ support_points[1] := [1,0]:
+ support_points[2] := [1,1]:
+ support_points[3] := [0,1]:
+ support_points[4] := [1/2,1/2]:
phi[0] := proc(x,y) if(y<1-x) then 1-x-y; else 0; fi; end:
phi[1] := proc(x,y) if(y<x) then x-y; else 0; fi; end:
- phi[2](x,y) - phi[3](x,y) ; end:
#points on children: let them be indexed one-based, as are
- #the ansatz_points
+ #the support_points
points[0] := array(0..n_functions-1, 1..2):
points[1] := array(0..n_functions-1, 1..2):
points[2] := array(0..n_functions-1, 1..2):
points[3] := array(0..n_functions-1, 1..2):
for i from 0 to n_functions-1 do
- points[0][i,1] := ansatz_points[i][1]/2:
- points[0][i,2] := ansatz_points[i][2]/2:
+ points[0][i,1] := support_points[i][1]/2:
+ points[0][i,2] := support_points[i][2]/2:
- points[1][i,1] := ansatz_points[i][1]/2+1/2:
- points[1][i,2] := ansatz_points[i][2]/2:
+ points[1][i,1] := support_points[i][1]/2+1/2:
+ points[1][i,2] := support_points[i][2]/2:
- points[2][i,1] := ansatz_points[i][1]/2+1/2:
- points[2][i,2] := ansatz_points[i][2]/2+1/2:
+ points[2][i,1] := support_points[i][1]/2+1/2:
+ points[2][i,2] := support_points[i][2]/2+1/2:
- points[3][i,1] := ansatz_points[i][1]/2:
- points[3][i,2] := ansatz_points[i][2]/2+1/2:
+ points[3][i,1] := support_points[i][1]/2:
+ points[3][i,2] := support_points[i][2]/2+1/2:
od:
prolongation := array(0..3,0..n_functions-1, 0..n_functions-1):
phi_eta[4] := proc(x,y) 1 - phi_eta[0](x,y) - phi_eta[1](x,y)
- phi_eta[2](x,y) - phi_eta[3](x,y) ; end:
- # define an array of the ansatz points in real space; the first
+ # define an array of the support points in real space; the first
# four are the vertices, the last one is the crossing point of
# the two diagonals
print ("Computing cross point"):
template <>
-void FECrissCross<1>::get_unit_ansatz_points (vector<Point<1> >&) const {
+void FECrissCross<1>::get_unit_support_points (vector<Point<1> >&) const {
Assert (false, ExcNotUseful());
};
template <>
-void FECrissCross<1>::get_ansatz_points (const DoFHandler<1>::cell_iterator &,
+void FECrissCross<1>::get_support_points (const DoFHandler<1>::cell_iterator &,
const Boundary<1> &,
vector<Point<1> > &) const {
Assert (false, ExcNotUseful());
template <>
-void FECrissCross<1>::get_face_ansatz_points (const DoFHandler<1>::face_iterator &,
+void FECrissCross<1>::get_face_support_points (const DoFHandler<1>::face_iterator &,
const Boundary<1> &,
vector<Point<1> > &) const {
Assert (false, ExcNotUseful());
template <>
-void FECrissCross<2>::get_unit_ansatz_points (vector<Point<2> > &unit_points) const {
+void FECrissCross<2>::get_unit_support_points (vector<Point<2> > &unit_points) const {
Assert(unit_points.size()==total_dofs,
ExcWrongFieldDimension (unit_points.size(), total_dofs));
template <>
-void FECrissCross<2>::get_ansatz_points (const DoFHandler<2>::cell_iterator &cell,
+void FECrissCross<2>::get_support_points (const DoFHandler<2>::cell_iterator &cell,
const Boundary<2> &,
- vector<Point<2> > &ansatz_points) const {
+ vector<Point<2> > &support_points) const {
const unsigned int dim = 2;
- Assert (ansatz_points.size() == total_dofs,
- ExcWrongFieldDimension (ansatz_points.size(), total_dofs));
+ Assert (support_points.size() == total_dofs,
+ ExcWrongFieldDimension (support_points.size(), total_dofs));
// copy vertices
for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell; ++vertex)
- ansatz_points[vertex] = cell->vertex(vertex);
+ support_points[vertex] = cell->vertex(vertex);
/*
- last ansatz point is the common point of the two diagonals; the formula for
+ last support point is the common point of the two diagonals; the formula for
the computation is a bit lengthy but straightforward. You can get it with
the small Maple script printed at the beginning of this file.
*/
const double t5 = x3*y0;
const double t14 = (t1-t2+x1*y3-t4+t5-x3*y1)/(t1-t2-x2*y1+x2*y3-t4+x1*y2+t5-x3*y2);
const double t15 = 1.0-t14;
- ansatz_points[4](0) = t15*x0+t14*x2;
- ansatz_points[4](1) = t15*y0+t14*y2;
+ support_points[4](0) = t15*x0+t14*x2;
+ support_points[4](1) = t15*y0+t14*y2;
};
template <>
-void FECrissCross<2>::get_face_ansatz_points (const DoFHandler<2>::face_iterator &face,
+void FECrissCross<2>::get_face_support_points (const DoFHandler<2>::face_iterator &face,
const Boundary<2> &,
- vector<Point<2> > &ansatz_points) const {
+ vector<Point<2> > &support_points) const {
const unsigned int dim = 2;
- Assert ((ansatz_points.size() == dofs_per_face) &&
- (ansatz_points.size() == GeometryInfo<dim>::vertices_per_face),
- ExcWrongFieldDimension (ansatz_points.size(),
+ Assert ((support_points.size() == dofs_per_face) &&
+ (support_points.size() == GeometryInfo<dim>::vertices_per_face),
+ ExcWrongFieldDimension (support_points.size(),
GeometryInfo<dim>::vertices_per_face));
for (unsigned int vertex=0; vertex<dofs_per_face; ++vertex)
- ansatz_points[vertex] = face->vertex(vertex);
+ support_points[vertex] = face->vertex(vertex);
};
const vector<Point<dim> > &unit_points,
vector<dFMatrix> &jacobians,
const bool compute_jacobians,
- vector<Point<dim> > &ansatz_points,
+ vector<Point<dim> > &support_points,
const bool,
vector<Point<dim> > &q_points,
const bool compute_q_points,
ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
Assert (q_points.size() == unit_points.size(),
ExcWrongFieldDimension(q_points.size(), unit_points.size()));
- Assert (ansatz_points.size() == total_dofs,
- ExcWrongFieldDimension(ansatz_points.size(), total_dofs));
+ Assert (support_points.size() == total_dofs,
+ ExcWrongFieldDimension(support_points.size(), total_dofs));
unsigned int n_points=unit_points.size();
- // we need the ansatz points in any
+ // we need the support points in any
// way, wanted or not by the user
- get_ansatz_points (cell, boundary, ansatz_points);
+ get_support_points (cell, boundary, support_points);
if (compute_q_points)
{
// x_l(xi_l) = sum_j p_j N_j(xi_l)
for (unsigned int j=0; j<n_transform_functions; ++j)
for (unsigned int l=0; l<n_points; ++l)
- q_points[l] += ansatz_points[j] * shape_values_transform(j, l);
+ q_points[l] += support_points[j] * shape_values_transform(j, l);
};
i=0..dim-1
j=0..dim-1
M_{ij}(l) = 0
- s=0..n_ansatz_points
+ s=0..n_support_points
M_{ij}(l) += p_i(s) d(N_s(l))/d(xi_j)
However, we rewrite the loops to only compute the gradient once for
const Point<dim> gradient = shape_grad_transform[s][l];
for (unsigned int i=0; i<dim; ++i)
for (unsigned int j=0; j<dim; ++j)
- M(i,j) += ansatz_points[s](i) * gradient(j);
+ M(i,j) += support_points[s](i) * gradient(j);
};
jacobians[l].invert(M);
};
--------------------------------------------------------------------------
n_functions := 4;
- ansatz_points := array(0..n_functions-1);
- ansatz_points[0] := 0;
- ansatz_points[1] := 1;
- ansatz_points[2] := 1/3;
- ansatz_points[3] := 2/3;
+ 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);
od;
values[i+1] := 1;
- shifted_ansatz_points := array (1..n_functions);
+ shifted_support_points := array (1..n_functions);
for j from 1 to n_functions do
- shifted_ansatz_points[j] := ansatz_points[j-1];
+ shifted_support_points[j] := support_points[j-1];
od;
- phi_polynom[i] := interp (shifted_ansatz_points, values, xi);
+ phi_polynom[i] := interp (shifted_support_points, values, xi);
grad_phi_polynom[i] := diff(phi_polynom[i], xi);
od;
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] := ansatz_points[i]/2;
- points[1][i] := ansatz_points[i]/2+1/2;
+ 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);
n_functions := 16:
n_face_functions := 4:
- ansatz_function := (a1 + a2*xi + a3*xi*xi + a4*xi*xi*xi) +
+ trial_function := (a1 + a2*xi + a3*xi*xi + a4*xi*xi*xi) +
(b1 + b2*xi + b3*xi*xi + b4*xi*xi*xi)*eta +
(c1 + c2*xi + c3*xi*xi + c4*xi*xi*xi)*eta*eta +
(d1 + d2*xi + d3*xi*xi + d4*xi*xi*xi)*eta*eta*eta:
- face_ansatz_function := a + b*xi + c*xi*xi + d*xi*xi*xi:
- # note: ansatz_points[i] is a vector which is indexed from
+ face_trial_function := a + b*xi + c*xi*xi + d*xi*xi*xi:
+ # note: support_points[i] is a vector which is indexed from
# one and not from zero!
- ansatz_points := array(0..n_functions-1):
- ansatz_points[0] := [0,0]:
- ansatz_points[1] := [1,0]:
- ansatz_points[2] := [1,1]:
- ansatz_points[3] := [0,1]:
- ansatz_points[4] := [1/3,0]:
- ansatz_points[5] := [2/3,0]:
- ansatz_points[6] := [1,1/3]:
- ansatz_points[7] := [1,2/3]:
- ansatz_points[8] := [1/3,1]:
- ansatz_points[9] := [2/3,1]:
- ansatz_points[10]:= [0,1/3]:
- ansatz_points[11]:= [0,2/3]:
- ansatz_points[12]:= [1/3,1/3]:
- ansatz_points[13]:= [2/3,1/3]:
- ansatz_points[14]:= [2/3,2/3]:
- ansatz_points[15]:= [1/3,2/3]:
-
- face_ansatz_points := array(0..n_face_functions-1):
- face_ansatz_points[0] := 0:
- face_ansatz_points[1] := 1:
- face_ansatz_points[2] := 1/3:
- face_ansatz_points[3] := 2/3:
- constrained_face_ansatz_points := array(0..2*(n_face_functions-2)+1-1):
- constrained_face_ansatz_points[0] := 1/2:
- constrained_face_ansatz_points[1] := 1/6:
- constrained_face_ansatz_points[2] := 2/6:
- constrained_face_ansatz_points[3] := 4/6:
- constrained_face_ansatz_points[4] := 5/6:
+ support_points := array(0..n_functions-1):
+ support_points[0] := [0,0]:
+ support_points[1] := [1,0]:
+ support_points[2] := [1,1]:
+ support_points[3] := [0,1]:
+ support_points[4] := [1/3,0]:
+ support_points[5] := [2/3,0]:
+ support_points[6] := [1,1/3]:
+ support_points[7] := [1,2/3]:
+ support_points[8] := [1/3,1]:
+ support_points[9] := [2/3,1]:
+ support_points[10]:= [0,1/3]:
+ support_points[11]:= [0,2/3]:
+ support_points[12]:= [1/3,1/3]:
+ support_points[13]:= [2/3,1/3]:
+ support_points[14]:= [2/3,2/3]:
+ support_points[15]:= [1/3,2/3]:
+
+ face_support_points := array(0..n_face_functions-1):
+ face_support_points[0] := 0:
+ face_support_points[1] := 1:
+ face_support_points[2] := 1/3:
+ face_support_points[3] := 2/3:
+ constrained_face_support_points := array(0..2*(n_face_functions-2)+1-1):
+ constrained_face_support_points[0] := 1/2:
+ constrained_face_support_points[1] := 1/6:
+ constrained_face_support_points[2] := 2/6:
+ constrained_face_support_points[3] := 4/6:
+ constrained_face_support_points[4] := 5/6:
phi_polynom := array(0..n_functions-1):
grad_phi_polynom := array(0..n_functions-1,0..1):
equation_system := {}:
for j from 0 to n_functions-1 do
- poly := subs(xi=ansatz_points[j][1],
- eta=ansatz_points[j][2],
- ansatz_function):
+ poly := subs(xi=support_points[j][1],
+ eta=support_points[j][2],
+ trial_function):
if (i=j) then
equation_system := equation_system union {poly = 1}:
else
fi:
od:
- phi_polynom[i] := subs(solve(equation_system), ansatz_function):
+ phi_polynom[i] := subs(solve(equation_system), trial_function):
grad_phi_polynom[i,0] := diff(phi_polynom[i], xi):
grad_phi_polynom[i,1] := diff(phi_polynom[i], eta):
od:
#points on children: let them be indexed one-based, as are
- #the ansatz_points
+ #the support_points
points[0] := array(0..n_functions-1, 1..2):
points[1] := array(0..n_functions-1, 1..2):
points[2] := array(0..n_functions-1, 1..2):
points[3] := array(0..n_functions-1, 1..2):
for i from 0 to n_functions-1 do
- points[0][i,1] := ansatz_points[i][1]/2:
- points[0][i,2] := ansatz_points[i][2]/2:
+ points[0][i,1] := support_points[i][1]/2:
+ points[0][i,2] := support_points[i][2]/2:
- points[1][i,1] := ansatz_points[i][1]/2+1/2:
- points[1][i,2] := ansatz_points[i][2]/2:
+ points[1][i,1] := support_points[i][1]/2+1/2:
+ points[1][i,2] := support_points[i][2]/2:
- points[2][i,1] := ansatz_points[i][1]/2+1/2:
- points[2][i,2] := ansatz_points[i][2]/2+1/2:
+ points[2][i,1] := support_points[i][1]/2+1/2:
+ points[2][i,2] := support_points[i][2]/2+1/2:
- points[3][i,1] := ansatz_points[i][1]/2:
- points[3][i,2] := ansatz_points[i][2]/2+1/2:
+ points[3][i,1] := support_points[i][1]/2:
+ points[3][i,2] := support_points[i][2]/2+1/2:
od:
print ("Computing prolongation matrices"):
od:
od:
- print ("computing ansatz points in real space"):
+ print ("computing support points in real space"):
for i from 0 to n_functions-1 do
- real_points[i,0] := subs(xi=ansatz_points[i][1],
- eta=ansatz_points[i][2], x_real);
- real_points[i,1] := subs(xi=ansatz_points[i][1],
- eta=ansatz_points[i][2], y_real);
+ real_points[i,0] := subs(xi=support_points[i][1],
+ eta=support_points[i][2], x_real);
+ real_points[i,1] := subs(xi=support_points[i][1],
+ eta=support_points[i][2], y_real);
od:
print ("computing interface constraint matrices"):
od:
values[i+1] := 1:
- shifted_face_ansatz_points := array (1..n_face_functions):
+ shifted_face_support_points := array (1..n_face_functions):
for j from 1 to n_face_functions do
- shifted_face_ansatz_points[j] := face_ansatz_points[j-1]:
+ shifted_face_support_points[j] := face_support_points[j-1]:
od:
- face_phi_polynom[i] := interp (shifted_face_ansatz_points, values, xi):
+ face_phi_polynom[i] := interp (shifted_face_support_points, values, xi):
od:
for i from 0 to 2*(n_face_functions-2)+1-1 do
for j from 0 to n_face_functions-1 do
- interface_constraints[i,j] := subs(xi=constrained_face_ansatz_points[i],
+ interface_constraints[i,j] := subs(xi=constrained_face_support_points[i],
face_phi_polynom[j]);
od:
od:
template <>
-void FECubicSub<1>::get_unit_ansatz_points (vector<Point<1> > &unit_points) const {
- FiniteElement<1>::get_unit_ansatz_points (unit_points);
+void FECubicSub<1>::get_unit_support_points (vector<Point<1> > &unit_points) const {
+ FiniteElement<1>::get_unit_support_points (unit_points);
};
template <>
-void FECubicSub<1>::get_ansatz_points (const typename DoFHandler<1>::cell_iterator &cell,
+void FECubicSub<1>::get_support_points (const typename DoFHandler<1>::cell_iterator &cell,
const Boundary<1> &boundary,
- vector<Point<1> > &ansatz_points) const {
- FiniteElement<1>::get_ansatz_points (cell, boundary, ansatz_points);
+ vector<Point<1> > &support_points) const {
+ FiniteElement<1>::get_support_points (cell, boundary, support_points);
};
template <>
-void FECubicSub<1>::get_face_ansatz_points (const typename DoFHandler<1>::face_iterator &,
+void FECubicSub<1>::get_face_support_points (const typename DoFHandler<1>::face_iterator &,
const Boundary<1> &,
vector<Point<1> > &) const {
Assert (false, ExcInternalError());
template <>
-void FECubicSub<2>::get_unit_ansatz_points (vector<Point<2> > &unit_points) const {
+void FECubicSub<2>::get_unit_support_points (vector<Point<2> > &unit_points) const {
Assert (unit_points.size() == total_dofs,
ExcWrongFieldDimension (unit_points.size(), total_dofs));
template <>
-void FECubicSub<2>::get_ansatz_points (const typename DoFHandler<2>::cell_iterator &cell,
+void FECubicSub<2>::get_support_points (const typename DoFHandler<2>::cell_iterator &cell,
const Boundary<2>&,
- vector<Point<2> > &ansatz_points) const {
- Assert (ansatz_points.size() == total_dofs,
- ExcWrongFieldDimension (ansatz_points.size(), total_dofs));
+ vector<Point<2> > &support_points) const {
+ Assert (support_points.size() == total_dofs,
+ ExcWrongFieldDimension (support_points.size(), total_dofs));
const double x[4] = { cell->vertex(0)(0),
cell->vertex(1)(0),
const double t45 = 2.0/9.0*x[2];
const double t48 = 2.0/9.0*y[0];
const double t50 = 2.0/9.0*y[2];
- ansatz_points[0](0) = x[0];
- ansatz_points[0](1) = y[0];
- ansatz_points[1](0) = x[1];
- ansatz_points[1](1) = y[1];
- ansatz_points[2](0) = x[2];
- ansatz_points[2](1) = y[2];
- ansatz_points[3](0) = x[3];
- ansatz_points[3](1) = y[3];
- ansatz_points[4](0) = t1+t2;
- ansatz_points[4](1) = t4+t5;
- ansatz_points[5](0) = t7+t8;
- ansatz_points[5](1) = t10+t11;
- ansatz_points[6](0) = t8+t13;
- ansatz_points[6](1) = t11+t15;
- ansatz_points[7](0) = t2+t17;
- ansatz_points[7](1) = t5+t19;
- ansatz_points[8](0) = t13+t21;
- ansatz_points[8](1) = t15+t23;
- ansatz_points[9](0) = t17+t25;
- ansatz_points[9](1) = t19+t27;
- ansatz_points[10](0) = t1+t25;
- ansatz_points[10](1) = t4+t27;
- ansatz_points[11](0) = t7+t21;
- ansatz_points[11](1) = t10+t23;
- ansatz_points[12](0) = 4.0/9.0*x[0]+t34+x[2]/9.0+t36;
- ansatz_points[12](1) = 4.0/9.0*y[0]+t39+y[2]/9.0+t41;
- ansatz_points[13](0) = t43+4.0/9.0*x[1]+t45+x[3]/9.0;
- ansatz_points[13](1) = t48+4.0/9.0*y[1]+t50+y[3]/9.0;
- ansatz_points[14](0) = x[0]/9.0+t34+4.0/9.0*x[2]+t36;
- ansatz_points[14](1) = y[0]/9.0+t39+4.0/9.0*y[2]+t41;
- ansatz_points[15](0) = t43+x[1]/9.0+t45+4.0/9.0*x[3];
- ansatz_points[15](1) = t48+y[1]/9.0+t50+4.0/9.0*y[3];
+ support_points[0](0) = x[0];
+ support_points[0](1) = y[0];
+ support_points[1](0) = x[1];
+ support_points[1](1) = y[1];
+ support_points[2](0) = x[2];
+ support_points[2](1) = y[2];
+ support_points[3](0) = x[3];
+ support_points[3](1) = y[3];
+ support_points[4](0) = t1+t2;
+ support_points[4](1) = t4+t5;
+ support_points[5](0) = t7+t8;
+ support_points[5](1) = t10+t11;
+ support_points[6](0) = t8+t13;
+ support_points[6](1) = t11+t15;
+ support_points[7](0) = t2+t17;
+ support_points[7](1) = t5+t19;
+ support_points[8](0) = t13+t21;
+ support_points[8](1) = t15+t23;
+ support_points[9](0) = t17+t25;
+ support_points[9](1) = t19+t27;
+ support_points[10](0) = t1+t25;
+ support_points[10](1) = t4+t27;
+ support_points[11](0) = t7+t21;
+ support_points[11](1) = t10+t23;
+ support_points[12](0) = 4.0/9.0*x[0]+t34+x[2]/9.0+t36;
+ support_points[12](1) = 4.0/9.0*y[0]+t39+y[2]/9.0+t41;
+ support_points[13](0) = t43+4.0/9.0*x[1]+t45+x[3]/9.0;
+ support_points[13](1) = t48+4.0/9.0*y[1]+t50+y[3]/9.0;
+ support_points[14](0) = x[0]/9.0+t34+4.0/9.0*x[2]+t36;
+ support_points[14](1) = y[0]/9.0+t39+4.0/9.0*y[2]+t41;
+ support_points[15](0) = t43+x[1]/9.0+t45+4.0/9.0*x[3];
+ support_points[15](1) = t48+y[1]/9.0+t50+4.0/9.0*y[3];
};
template <>
-void FECubicSub<2>::get_face_ansatz_points (const typename DoFHandler<2>::face_iterator &face,
+void FECubicSub<2>::get_face_support_points (const typename DoFHandler<2>::face_iterator &face,
const Boundary<2> &,
- vector<Point<2> > &ansatz_points) const {
- Assert (ansatz_points.size() == dofs_per_face,
- ExcWrongFieldDimension (ansatz_points.size(), dofs_per_face));
+ vector<Point<2> > &support_points) const {
+ Assert (support_points.size() == dofs_per_face,
+ ExcWrongFieldDimension (support_points.size(), dofs_per_face));
for (unsigned int vertex=0; vertex<2; ++vertex)
- ansatz_points[vertex] = face->vertex(vertex);
- ansatz_points[2] = (2*ansatz_points[0] + ansatz_points[1]) / 3;
- ansatz_points[3] = (ansatz_points[0] + 2*ansatz_points[1]) / 3;
+ support_points[vertex] = face->vertex(vertex);
+ support_points[2] = (2*support_points[0] + support_points[1]) / 3;
+ support_points[3] = (support_points[0] + 2*support_points[1]) / 3;
};
template <>
-void FELinear<1>::get_unit_ansatz_points (vector<Point<1> > &ansatz_points) const {
- FiniteElement<1>::get_unit_ansatz_points (ansatz_points);
+void FELinear<1>::get_unit_support_points (vector<Point<1> > &support_points) const {
+ FiniteElement<1>::get_unit_support_points (support_points);
};
template <>
-void FELinear<1>::get_ansatz_points (const typename DoFHandler<1>::cell_iterator &cell,
+void FELinear<1>::get_support_points (const typename DoFHandler<1>::cell_iterator &cell,
const Boundary<1> &boundary,
- vector<Point<1> > &ansatz_points) const {
- FiniteElement<1>::get_ansatz_points (cell, boundary, ansatz_points);
+ vector<Point<1> > &support_points) const {
+ FiniteElement<1>::get_support_points (cell, boundary, support_points);
};
template <>
-void FELinear<1>::get_face_ansatz_points (const typename DoFHandler<1>::face_iterator &,
+void FELinear<1>::get_face_support_points (const typename DoFHandler<1>::face_iterator &,
const Boundary<1> &,
vector<Point<1> > &) const {
Assert (false, ExcInternalError());
template <>
-void FELinear<2>::get_unit_ansatz_points (vector<Point<2> > &unit_points) const {
+void FELinear<2>::get_unit_support_points (vector<Point<2> > &unit_points) const {
Assert (unit_points.size() == total_dofs,
ExcWrongFieldDimension (unit_points.size(), total_dofs));
template <int dim>
-void FELinear<dim>::get_ansatz_points (const typename DoFHandler<dim>::cell_iterator &cell,
+void FELinear<dim>::get_support_points (const typename DoFHandler<dim>::cell_iterator &cell,
const Boundary<dim> &,
- vector<Point<dim> > &ansatz_points) const {
- Assert (ansatz_points.size() == total_dofs,
- ExcWrongFieldDimension (ansatz_points.size(), total_dofs));
+ vector<Point<dim> > &support_points) const {
+ Assert (support_points.size() == total_dofs,
+ ExcWrongFieldDimension (support_points.size(), total_dofs));
for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell; ++vertex)
- ansatz_points[vertex] = cell->vertex(vertex);
+ support_points[vertex] = cell->vertex(vertex);
};
template <int dim>
-void FELinear<dim>::get_face_ansatz_points (const typename DoFHandler<dim>::face_iterator &face,
+void FELinear<dim>::get_face_support_points (const typename DoFHandler<dim>::face_iterator &face,
const Boundary<dim> &,
- vector<Point<dim> > &ansatz_points) const {
- Assert ((ansatz_points.size() == dofs_per_face) &&
- (ansatz_points.size() == GeometryInfo<dim>::vertices_per_face),
- ExcWrongFieldDimension (ansatz_points.size(),
+ vector<Point<dim> > &support_points) const {
+ Assert ((support_points.size() == dofs_per_face) &&
+ (support_points.size() == GeometryInfo<dim>::vertices_per_face),
+ ExcWrongFieldDimension (support_points.size(),
GeometryInfo<dim>::vertices_per_face));
for (unsigned int vertex=0; vertex<dofs_per_face; ++vertex)
- ansatz_points[vertex] = face->vertex(vertex);
+ support_points[vertex] = face->vertex(vertex);
};
template <>
-void FEQuadraticSub<1>::get_unit_ansatz_points (vector<Point<1> > &unit_points) const {
- FiniteElement<1>::get_unit_ansatz_points (unit_points);
+void FEQuadraticSub<1>::get_unit_support_points (vector<Point<1> > &unit_points) const {
+ FiniteElement<1>::get_unit_support_points (unit_points);
};
template <>
-void FEQuadraticSub<1>::get_ansatz_points (const typename DoFHandler<1>::cell_iterator &cell,
+void FEQuadraticSub<1>::get_support_points (const typename DoFHandler<1>::cell_iterator &cell,
const Boundary<1> &boundary,
- vector<Point<1> > &ansatz_points) const {
- FiniteElement<1>::get_ansatz_points (cell, boundary, ansatz_points);
+ vector<Point<1> > &support_points) const {
+ FiniteElement<1>::get_support_points (cell, boundary, support_points);
};
template <>
-void FEQuadraticSub<1>::get_face_ansatz_points (const typename DoFHandler<1>::face_iterator &,
+void FEQuadraticSub<1>::get_face_support_points (const typename DoFHandler<1>::face_iterator &,
const Boundary<1> &,
vector<Point<1> > &) const {
Assert (false, ExcInternalError());
template <>
-void FEQuadraticSub<2>::get_unit_ansatz_points (vector<Point<2> > &unit_points) const {
+void FEQuadraticSub<2>::get_unit_support_points (vector<Point<2> > &unit_points) const {
Assert (unit_points.size() == total_dofs,
ExcWrongFieldDimension (unit_points.size(), total_dofs));
template <>
-void FEQuadraticSub<2>::get_ansatz_points (const typename DoFHandler<2>::cell_iterator &cell,
+void FEQuadraticSub<2>::get_support_points (const typename DoFHandler<2>::cell_iterator &cell,
const Boundary<2>&,
- vector<Point<2> > &ansatz_points) const {
- Assert (ansatz_points.size() == total_dofs,
- ExcWrongFieldDimension (ansatz_points.size(), total_dofs));
+ vector<Point<2> > &support_points) const {
+ Assert (support_points.size() == total_dofs,
+ ExcWrongFieldDimension (support_points.size(), total_dofs));
for (unsigned int vertex=0; vertex<4; ++vertex)
- ansatz_points[vertex] = cell->vertex(vertex);
+ support_points[vertex] = cell->vertex(vertex);
// for the bilinear mapping, the centers
// of the face on the unit cell are mapped
// to the mean coordinates of the vertices
for (unsigned int line=0; line<4; ++line)
- ansatz_points[4+line] = (cell->line(line)->vertex(0) +
+ support_points[4+line] = (cell->line(line)->vertex(0) +
cell->line(line)->vertex(1)) / 2;
// same for the center of the square:
// since all four linear basis functions
// take on the value 1/4 at the center,
// the center is mapped to the mean
// coordinates of the four vertices
- ansatz_points[8] = (ansatz_points[0] +
- ansatz_points[1] +
- ansatz_points[2] +
- ansatz_points[3]) / 4;
+ support_points[8] = (support_points[0] +
+ support_points[1] +
+ support_points[2] +
+ support_points[3]) / 4;
};
template <>
-void FEQuadraticSub<2>::get_face_ansatz_points (const typename DoFHandler<2>::face_iterator &face,
+void FEQuadraticSub<2>::get_face_support_points (const typename DoFHandler<2>::face_iterator &face,
const Boundary<2> &,
- vector<Point<2> > &ansatz_points) const {
- Assert (ansatz_points.size() == dofs_per_face,
- ExcWrongFieldDimension (ansatz_points.size(), dofs_per_face));
+ vector<Point<2> > &support_points) const {
+ Assert (support_points.size() == dofs_per_face,
+ ExcWrongFieldDimension (support_points.size(), dofs_per_face));
for (unsigned int vertex=0; vertex<2; ++vertex)
- ansatz_points[vertex] = face->vertex(vertex);
- ansatz_points[2] = (ansatz_points[0] + ansatz_points[1]) / 2;
+ support_points[vertex] = face->vertex(vertex);
+ support_points[2] = (support_points[0] + support_points[1]) / 2;
};
--------------------------------------------------------------------------
n_functions := 5;
- ansatz_points := array(0..n_functions-1);
- ansatz_points[0] := 0;
- ansatz_points[1] := 1;
- ansatz_points[2] := 1/4;
- ansatz_points[3] := 2/4;
- ansatz_points[4] := 3/4;
+ 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);
od;
values[i+1] := 1;
- shifted_ansatz_points := array (1..n_functions);
+ shifted_support_points := array (1..n_functions);
for j from 1 to n_functions do
- shifted_ansatz_points[j] := ansatz_points[j-1];
+ shifted_support_points[j] := support_points[j-1];
od;
- phi_polynom[i] := interp (shifted_ansatz_points, values, xi);
+ phi_polynom[i] := interp (shifted_support_points, values, xi);
grad_phi_polynom[i] := diff(phi_polynom[i], xi);
od;
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] := ansatz_points[i]/2;
- points[1][i] := ansatz_points[i]/2+1/2;
+ 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);
n_functions := 25:
n_face_functions := 5:
- ansatz_function := (a1 + a2*xi + a3*xi*xi + a4*xi**3 + a5*xi**4) +
+ trial_function := (a1 + a2*xi + a3*xi*xi + a4*xi**3 + a5*xi**4) +
(b1 + b2*xi + b3*xi*xi + b4*xi**3 + b5*xi**4)*eta +
(c1 + c2*xi + c3*xi*xi + c4*xi**3 + c5*xi**4)*eta*eta +
(d1 + d2*xi + d3*xi*xi + d4*xi**3 + d5*xi**4)*eta**3 +
(e1 + e2*xi + e3*xi*xi + e4*xi**3 + e5*xi**4)*eta**4:
- face_ansatz_function := a + b*xi + c*xi*xi + d*xi**3 + e*xi**4:
- # note: ansatz_points[i] is a vector which is indexed from
+ face_trial_function := a + b*xi + c*xi*xi + d*xi**3 + e*xi**4:
+ # note: support_points[i] is a vector which is indexed from
# one and not from zero!
- ansatz_points := array(0..n_functions-1):
- ansatz_points[0] := [0,0]:
- ansatz_points[1] := [1,0]:
- ansatz_points[2] := [1,1]:
- ansatz_points[3] := [0,1]:
- ansatz_points[4] := [1/4,0]:
- ansatz_points[5] := [2/4,0]:
- ansatz_points[6] := [3/4,0]:
- ansatz_points[7] := [1,1/4]:
- ansatz_points[8] := [1,2/4]:
- ansatz_points[9] := [1,3/4]:
- ansatz_points[10] := [1/4,1]:
- ansatz_points[11] := [2/4,1]:
- ansatz_points[12] := [3/4,1]:
- ansatz_points[13] := [0,1/4]:
- ansatz_points[14] := [0,2/4]:
- ansatz_points[15] := [0,3/4]:
- ansatz_points[16] := [1/4,1/4]:
- ansatz_points[17] := [3/4,1/4]:
- ansatz_points[18] := [3/4,3/4]:
- ansatz_points[19] := [1/4,3/4]:
- ansatz_points[20] := [1/2,1/4]:
- ansatz_points[21] := [3/4,1/2]:
- ansatz_points[22] := [1/2,3/4]:
- ansatz_points[23] := [1/4,1/2]:
- ansatz_points[24] := [1/2,1/2]:
-
- face_ansatz_points := array(0..n_face_functions-1):
- face_ansatz_points[0] := 0:
- face_ansatz_points[1] := 1:
- face_ansatz_points[2] := 1/4:
- face_ansatz_points[3] := 2/4:
- face_ansatz_points[4] := 3/4:
- constrained_face_ansatz_points := array(0..2*(n_face_functions-2)+1-1):
- constrained_face_ansatz_points[0] := 1/2:
- constrained_face_ansatz_points[1] := 1/8:
- constrained_face_ansatz_points[2] := 2/8:
- constrained_face_ansatz_points[3] := 3/8:
- constrained_face_ansatz_points[4] := 5/8:
- constrained_face_ansatz_points[5] := 6/8:
- constrained_face_ansatz_points[6] := 7/8:
+ support_points := array(0..n_functions-1):
+ support_points[0] := [0,0]:
+ support_points[1] := [1,0]:
+ support_points[2] := [1,1]:
+ support_points[3] := [0,1]:
+ support_points[4] := [1/4,0]:
+ support_points[5] := [2/4,0]:
+ support_points[6] := [3/4,0]:
+ support_points[7] := [1,1/4]:
+ support_points[8] := [1,2/4]:
+ support_points[9] := [1,3/4]:
+ support_points[10] := [1/4,1]:
+ support_points[11] := [2/4,1]:
+ support_points[12] := [3/4,1]:
+ support_points[13] := [0,1/4]:
+ support_points[14] := [0,2/4]:
+ support_points[15] := [0,3/4]:
+ support_points[16] := [1/4,1/4]:
+ support_points[17] := [3/4,1/4]:
+ support_points[18] := [3/4,3/4]:
+ support_points[19] := [1/4,3/4]:
+ support_points[20] := [1/2,1/4]:
+ support_points[21] := [3/4,1/2]:
+ support_points[22] := [1/2,3/4]:
+ support_points[23] := [1/4,1/2]:
+ support_points[24] := [1/2,1/2]:
+
+ face_support_points := array(0..n_face_functions-1):
+ face_support_points[0] := 0:
+ face_support_points[1] := 1:
+ face_support_points[2] := 1/4:
+ face_support_points[3] := 2/4:
+ face_support_points[4] := 3/4:
+ constrained_face_support_points := array(0..2*(n_face_functions-2)+1-1):
+ constrained_face_support_points[0] := 1/2:
+ constrained_face_support_points[1] := 1/8:
+ constrained_face_support_points[2] := 2/8:
+ constrained_face_support_points[3] := 3/8:
+ constrained_face_support_points[4] := 5/8:
+ constrained_face_support_points[5] := 6/8:
+ constrained_face_support_points[6] := 7/8:
phi_polynom := array(0..n_functions-1):
grad_phi_polynom := array(0..n_functions-1,0..1):
equation_system := {}:
for j from 0 to n_functions-1 do
- poly := subs(xi=ansatz_points[j][1],
- eta=ansatz_points[j][2],
- ansatz_function):
+ poly := subs(xi=support_points[j][1],
+ eta=support_points[j][2],
+ trial_function):
if (i=j) then
equation_system := equation_system union {poly = 1}:
else
fi:
od:
- phi_polynom[i] := subs(solve(equation_system), ansatz_function):
+ phi_polynom[i] := subs(solve(equation_system), trial_function):
grad_phi_polynom[i,0] := diff(phi_polynom[i], xi):
grad_phi_polynom[i,1] := diff(phi_polynom[i], eta):
od:
#points on children: let them be indexed one-based, as are
- #the ansatz_points
+ #the support_points
points[0] := array(0..n_functions-1, 1..2):
points[1] := array(0..n_functions-1, 1..2):
points[2] := array(0..n_functions-1, 1..2):
points[3] := array(0..n_functions-1, 1..2):
for i from 0 to n_functions-1 do
- points[0][i,1] := ansatz_points[i][1]/2:
- points[0][i,2] := ansatz_points[i][2]/2:
+ points[0][i,1] := support_points[i][1]/2:
+ points[0][i,2] := support_points[i][2]/2:
- points[1][i,1] := ansatz_points[i][1]/2+1/2:
- points[1][i,2] := ansatz_points[i][2]/2:
+ points[1][i,1] := support_points[i][1]/2+1/2:
+ points[1][i,2] := support_points[i][2]/2:
- points[2][i,1] := ansatz_points[i][1]/2+1/2:
- points[2][i,2] := ansatz_points[i][2]/2+1/2:
+ points[2][i,1] := support_points[i][1]/2+1/2:
+ points[2][i,2] := support_points[i][2]/2+1/2:
- points[3][i,1] := ansatz_points[i][1]/2:
- points[3][i,2] := ansatz_points[i][2]/2+1/2:
+ points[3][i,1] := support_points[i][1]/2:
+ points[3][i,2] := support_points[i][2]/2+1/2:
od:
print ("Computing prolongation matrices"):
od:
od:
- print ("computing ansatz points in real space"):
+ print ("computing support points in real space"):
for i from 0 to n_functions-1 do
- real_points[i,0] := subs(xi=ansatz_points[i][1],
- eta=ansatz_points[i][2], x_real);
- real_points[i,1] := subs(xi=ansatz_points[i][1],
- eta=ansatz_points[i][2], y_real);
+ real_points[i,0] := subs(xi=support_points[i][1],
+ eta=support_points[i][2], x_real);
+ real_points[i,1] := subs(xi=support_points[i][1],
+ eta=support_points[i][2], y_real);
od:
print ("computing interface constraint matrices"):
od:
values[i+1] := 1:
- shifted_face_ansatz_points := array (1..n_face_functions):
+ shifted_face_support_points := array (1..n_face_functions):
for j from 1 to n_face_functions do
- shifted_face_ansatz_points[j] := face_ansatz_points[j-1]:
+ shifted_face_support_points[j] := face_support_points[j-1]:
od:
- face_phi_polynom[i] := interp (shifted_face_ansatz_points, values, xi):
+ face_phi_polynom[i] := interp (shifted_face_support_points, values, xi):
od:
for i from 0 to 2*(n_face_functions-2)+1-1 do
for j from 0 to n_face_functions-1 do
- interface_constraints[i,j] := subs(xi=constrained_face_ansatz_points[i],
+ interface_constraints[i,j] := subs(xi=constrained_face_support_points[i],
face_phi_polynom[j]);
od:
od:
template <>
-void FEQuarticSub<1>::get_unit_ansatz_points (vector<Point<1> > &unit_points) const {
- FiniteElement<1>::get_unit_ansatz_points (unit_points);
+void FEQuarticSub<1>::get_unit_support_points (vector<Point<1> > &unit_points) const {
+ FiniteElement<1>::get_unit_support_points (unit_points);
};
template <>
-void FEQuarticSub<1>::get_ansatz_points (const typename DoFHandler<1>::cell_iterator &cell,
+void FEQuarticSub<1>::get_support_points (const typename DoFHandler<1>::cell_iterator &cell,
const Boundary<1> &boundary,
- vector<Point<1> > &ansatz_points) const {
- FiniteElement<1>::get_ansatz_points (cell, boundary, ansatz_points);
+ vector<Point<1> > &support_points) const {
+ FiniteElement<1>::get_support_points (cell, boundary, support_points);
};
template <>
-void FEQuarticSub<1>::get_face_ansatz_points (const typename DoFHandler<1>::face_iterator &,
+void FEQuarticSub<1>::get_face_support_points (const typename DoFHandler<1>::face_iterator &,
const Boundary<1> &,
vector<Point<1> > &) const {
Assert (false, ExcInternalError());
template <>
-void FEQuarticSub<2>::get_unit_ansatz_points (vector<Point<2> > &unit_points) const {
+void FEQuarticSub<2>::get_unit_support_points (vector<Point<2> > &unit_points) const {
Assert (unit_points.size() == total_dofs,
ExcWrongFieldDimension (unit_points.size(), total_dofs));
template <>
-void FEQuarticSub<2>::get_ansatz_points (const typename DoFHandler<2>::cell_iterator &cell,
+void FEQuarticSub<2>::get_support_points (const typename DoFHandler<2>::cell_iterator &cell,
const Boundary<2>&,
- vector<Point<2> > &ansatz_points) const {
- Assert (ansatz_points.size() == total_dofs,
- ExcWrongFieldDimension (ansatz_points.size(), total_dofs));
+ vector<Point<2> > &support_points) const {
+ Assert (support_points.size() == total_dofs,
+ ExcWrongFieldDimension (support_points.size(), total_dofs));
const double x[4] = { cell->vertex(0)(0),
cell->vertex(1)(0),
const double t90 = 3.0/8.0*x[3];
const double t92 = y[1]/8.0;
const double t93 = 3.0/8.0*y[3];
- ansatz_points[0](0) = x[0];
- ansatz_points[0](1) = y[0];
- ansatz_points[1](0) = x[1];
- ansatz_points[1](1) = y[1];
- ansatz_points[2](0) = x[2];
- ansatz_points[2](1) = y[2];
- ansatz_points[3](0) = x[3];
- ansatz_points[3](1) = y[3];
- ansatz_points[4](0) = t1+t2;
- ansatz_points[4](1) = t4+t5;
- ansatz_points[5](0) = x[0]/2.0+x[1]/2.0;
- ansatz_points[5](1) = y[0]/2.0+y[1]/2.0;
- ansatz_points[6](0) = t9+t10;
- ansatz_points[6](1) = t12+t13;
- ansatz_points[7](0) = t10+t15;
- ansatz_points[7](1) = t13+t17;
- ansatz_points[8](0) = x[1]/2.0+x[2]/2.0;
- ansatz_points[8](1) = y[1]/2.0+y[2]/2.0;
- ansatz_points[9](0) = t2+t21;
- ansatz_points[9](1) = t5+t23;
- ansatz_points[10](0) = t15+t25;
- ansatz_points[10](1) = t17+t27;
- ansatz_points[11](0) = x[2]/2.0+x[3]/2.0;
- ansatz_points[11](1) = y[2]/2.0+y[3]/2.0;
- ansatz_points[12](0) = t21+t31;
- ansatz_points[12](1) = t23+t33;
- ansatz_points[13](0) = t1+t31;
- ansatz_points[13](1) = t4+t33;
- ansatz_points[14](0) = x[0]/2.0+x[3]/2.0;
- ansatz_points[14](1) = y[0]/2.0+y[3]/2.0;
- ansatz_points[15](0) = t9+t25;
- ansatz_points[15](1) = t12+t27;
- ansatz_points[16](0) = 9.0/16.0*x[0]+t42+x[2]/16.0+t44;
- ansatz_points[16](1) = 9.0/16.0*y[0]+t47+y[2]/16.0+t49;
- ansatz_points[17](0) = t51+9.0/16.0*x[1]+t53+x[3]/16.0;
- ansatz_points[17](1) = t56+9.0/16.0*y[1]+t58+y[3]/16.0;
- ansatz_points[18](0) = x[0]/16.0+t42+9.0/16.0*x[2]+t44;
- ansatz_points[18](1) = y[0]/16.0+t47+9.0/16.0*y[2]+t49;
- ansatz_points[19](0) = t51+x[1]/16.0+t53+9.0/16.0*x[3];
- ansatz_points[19](1) = t56+y[1]/16.0+t58+9.0/16.0*y[3];
- ansatz_points[20](0) = t73+t74+t75+t76;
- ansatz_points[20](1) = t78+t79+t80+t81;
- ansatz_points[21](0) = t83+t74+t84+t76;
- ansatz_points[21](1) = t86+t79+t87+t81;
- ansatz_points[22](0) = t83+t89+t84+t90;
- ansatz_points[22](1) = t86+t92+t87+t93;
- ansatz_points[23](0) = t73+t89+t75+t90;
- ansatz_points[23](1) = t78+t92+t80+t93;
- ansatz_points[24](0) = x[0]/4.0+x[1]/4.0+x[2]/4.0+x[3]/4.0;
- ansatz_points[24](1) = y[0]/4.0+y[1]/4.0+y[2]/4.0+y[3]/4.0;
+ support_points[0](0) = x[0];
+ support_points[0](1) = y[0];
+ support_points[1](0) = x[1];
+ support_points[1](1) = y[1];
+ support_points[2](0) = x[2];
+ support_points[2](1) = y[2];
+ support_points[3](0) = x[3];
+ support_points[3](1) = y[3];
+ support_points[4](0) = t1+t2;
+ support_points[4](1) = t4+t5;
+ support_points[5](0) = x[0]/2.0+x[1]/2.0;
+ support_points[5](1) = y[0]/2.0+y[1]/2.0;
+ support_points[6](0) = t9+t10;
+ support_points[6](1) = t12+t13;
+ support_points[7](0) = t10+t15;
+ support_points[7](1) = t13+t17;
+ support_points[8](0) = x[1]/2.0+x[2]/2.0;
+ support_points[8](1) = y[1]/2.0+y[2]/2.0;
+ support_points[9](0) = t2+t21;
+ support_points[9](1) = t5+t23;
+ support_points[10](0) = t15+t25;
+ support_points[10](1) = t17+t27;
+ support_points[11](0) = x[2]/2.0+x[3]/2.0;
+ support_points[11](1) = y[2]/2.0+y[3]/2.0;
+ support_points[12](0) = t21+t31;
+ support_points[12](1) = t23+t33;
+ support_points[13](0) = t1+t31;
+ support_points[13](1) = t4+t33;
+ support_points[14](0) = x[0]/2.0+x[3]/2.0;
+ support_points[14](1) = y[0]/2.0+y[3]/2.0;
+ support_points[15](0) = t9+t25;
+ support_points[15](1) = t12+t27;
+ support_points[16](0) = 9.0/16.0*x[0]+t42+x[2]/16.0+t44;
+ support_points[16](1) = 9.0/16.0*y[0]+t47+y[2]/16.0+t49;
+ support_points[17](0) = t51+9.0/16.0*x[1]+t53+x[3]/16.0;
+ support_points[17](1) = t56+9.0/16.0*y[1]+t58+y[3]/16.0;
+ support_points[18](0) = x[0]/16.0+t42+9.0/16.0*x[2]+t44;
+ support_points[18](1) = y[0]/16.0+t47+9.0/16.0*y[2]+t49;
+ support_points[19](0) = t51+x[1]/16.0+t53+9.0/16.0*x[3];
+ support_points[19](1) = t56+y[1]/16.0+t58+9.0/16.0*y[3];
+ support_points[20](0) = t73+t74+t75+t76;
+ support_points[20](1) = t78+t79+t80+t81;
+ support_points[21](0) = t83+t74+t84+t76;
+ support_points[21](1) = t86+t79+t87+t81;
+ support_points[22](0) = t83+t89+t84+t90;
+ support_points[22](1) = t86+t92+t87+t93;
+ support_points[23](0) = t73+t89+t75+t90;
+ support_points[23](1) = t78+t92+t80+t93;
+ support_points[24](0) = x[0]/4.0+x[1]/4.0+x[2]/4.0+x[3]/4.0;
+ support_points[24](1) = y[0]/4.0+y[1]/4.0+y[2]/4.0+y[3]/4.0;
};
template <>
-void FEQuarticSub<2>::get_face_ansatz_points (const typename DoFHandler<2>::face_iterator &face,
+void FEQuarticSub<2>::get_face_support_points (const typename DoFHandler<2>::face_iterator &face,
const Boundary<2> &,
- vector<Point<2> > &ansatz_points) const {
- Assert (ansatz_points.size() == dofs_per_face,
- ExcWrongFieldDimension (ansatz_points.size(), dofs_per_face));
+ vector<Point<2> > &support_points) const {
+ Assert (support_points.size() == dofs_per_face,
+ ExcWrongFieldDimension (support_points.size(), dofs_per_face));
for (unsigned int vertex=0; vertex<2; ++vertex)
- ansatz_points[vertex] = face->vertex(vertex);
- ansatz_points[2] = (3*ansatz_points[0] + ansatz_points[1]) / 4;
- ansatz_points[3] = (ansatz_points[0] + ansatz_points[1]) / 2;
- ansatz_points[4] = (ansatz_points[0] + 3*ansatz_points[1]) / 4;
+ support_points[vertex] = face->vertex(vertex);
+ support_points[2] = (3*support_points[0] + support_points[1]) / 4;
+ support_points[3] = (support_points[0] + support_points[1]) / 2;
+ support_points[4] = (support_points[0] + 3*support_points[1]) / 4;
};