\section{Parametric elements}
\label{s:parametric}

With parametric elements, problems can be solved on meshes which dimensions are not necessarily equal to the world dimension. Therefore, problems on arbitrary manifolds can be solved. Furthermore, the vertex coordinates of the mesh can be flexible. Hence, moving meshes can be implemented. 

In this section, we solve equation (\ref{eq:ellipt}) with $f=2x_0$ ($x_0$ is the first component of $x$) on a torus. Then we rotate the torus about the $y$-axis and solve the problem again.

The torus can be created by revolving a circle about an axis coplanar with the circle, which does not touch the circle. We call $r_1$ the radius of the revolved circle and $r_2$ the radius of the revolution, which is the distance of the center of the tube to the center of the torus. In Figure \ref{f:torus}, a torus with its two radii $r_1$ and $r_2$ is shown. 
\begin{figure}
\center
\includegraphics[width=0.7\textwidth]{fig/torus.pdf}
\caption{A torus and a halfed torus with the two radiuses $r_1$ and $r_2$.}
\label{f:torus}
\end{figure}

We create a torus with center $(0; 0; 0)$ and the rotation axis in $z$-direction $(0; 0; 1)$. The projection of a point $x_0$ on the torus is implemented by the following steps:
\begin{enumerate}
\item $x_1$ is the projection of $x_0$ on the $xy$-plane
\item $x_2 = x_1 \frac{r_1}{||x_1||}$. Projection of $x_1$ on the sphere with radius $r_1$ with center $0$. Thereby, $x_2$ is used as the center of a sphere with radius $r_2$.  
\item $x_3 = x_0 - x_2$. Move coordinate system into the center of the sphere with center $x_2$. Thereby, $x_3$ contains the coordinates of $x_0$ in this new coordinate system.
\item $x_4 = x_3 \frac{r_2}{||x_3||}$. Thereby, $x_4$ is the projection of $x_3$ on the sphere with radius $r_2$.
\item $x_5 = x_4 + x_2$. Thereby, $x_5$ contains the coordinates of $x_4$ in the original coordinate system. It is the projection of $x_0$ on the torus.
\end{enumerate}

\begin{table}
\center
\begin{tabular}{|cc|c|c|c|}
\hline 
& & {\bf 1d} & {\bf 2d} & {\bf 3d} \\
\hline 
{\bf source code} & \tt src/ & \multicolumn{3}{|c|}{\tt torus.cc} \\
\hline 
{\bf parameter file} & \tt init/ & \tt - & \tt - & \tt torus.dat.3d \\
\hline 
{\bf macro file} & \tt macro/ & \tt -  & \tt - & \tt torus.macro.3d \\
\hline 
{\bf output files} & \tt output/ & \multicolumn{3}{|c|}{\tt torus.mesh/.dat, rotation1.mesh/.dat, rotation2.mesh/dat} \\
\hline 
\end{tabular}
\caption{Files of the {\tt torus} example.}
\label{t:torus files}
\end{table}

\subsection{Source code}

First, we define the rotation about the $y$-axis, which is used later to rotate the whole torus and the right hand side function.

\begin{lstlisting}{}
class YRotation
{
public:
  static WorldVector<double>& rotate(WorldVector<double> &x, double angle) 
  {
    double x0 = x[0] * cos(angle) + x[2] * sin(angle);
    x[2] = -x[0] * sin(angle) + x[2] * cos(angle);
    x[0] = x0;
    return x;
  }
};
\end{lstlisting}

The right hand side function $f$ has to follow the rotation of the torus.

\begin{lstlisting}{}
class F : public AbstractFunction<double, WorldVector<double> >
{
public:
  MEMORY_MANAGED(F);

  F(int degree) 
    : AbstractFunction<double, WorldVector<double> >(degree),
      rotation(0.0)
  {};

  double operator()(const WorldVector<double>& x) const {
    WorldVector<double> myX = x;
    YRotation::rotate(myX, -rotation);
    return = -2.0 * myX[0];
  }

  void rotate(double r) { rotation += r; };

private:
  double rotation;
};
\end{lstlisting}

Every time, the mesh is rotated, the right hand side function will be informed over the method \verb+rotate+. 

Now, we implement the projection on the torus.

\begin{lstlisting}{}
class TorusProject : public Projection
{
public:
  TorusProject(int id, 
	       ProjectionType type,
	       double radius1,
	       double radius2) 
    : Projection(id, type),
      radius1_(radius1),
      radius2_(radius2)
  {};

  virtual ~TorusProject() {};

  void project(WorldVector<double> &x) {

    WorldVector<double> xPlane = x;
    xPlane[2] = 0.0;

    double norm = std::sqrt(xPlane*xPlane);
    TEST_EXIT(norm != 0.0)("can't project vector x\n");
    
    WorldVector<double> center = xPlane;
    center *= radius1_ / norm;

    x -= center;

    norm = std::sqrt(x*x);
    TEST_EXIT(norm != 0.0)("can't project vector x\n");
    x *= radius2_/norm;

    x += center;
  };

protected:
  double radius1_;
  double radius2_;
};
\end{lstlisting}

In the main program, we create a torus projection as \verb+VOLUME_PROJECTION+ with ID 1. The values of $r_1$ and $r_2$ are chosen, such that the resulting torus is completely surrounded by the macro mesh that is defined later. 

\begin{lstlisting}{}
int main(int argc, char* argv[])
{
  FUNCNAME("torus main");

  // ===== check for init file =====
  TEST_EXIT(argc == 2)("usage: torus initfile\n");

  // ===== init parameters =====
  Parameters::init(false, argv[1]);

  // ===== create projection =====
  double r2 = (1.5 - 1.0/std::sqrt(2.0)) / 2.0;
  double r1 = 1.0/std::sqrt(2.0) + r2;

  NEW TorusProject(1, 
		   VOLUME_PROJECTION, 
		   r1,
		   r2);

  ...

  adapt->adapt();

  torus.writeFiles(adaptInfo, true);
\end{lstlisting}

The problem definition and the creation of the adaptation loop are done in the usual way (here, replaced by \verb+...+) . After the adaptation loop has returned, we write the result. 

Before we let the torus rotate, some variables are defined. We set the rotation angle to $\frac{\Pi}{3}$.

\begin{lstlisting}{}
  double rotation = M_PI/3.0;
  int i, j;
  int dim = torus.getMesh()->getDim();
  int dow = Global::getGeo(WORLD);

  DegreeOfFreedom dof;
  WorldVector<double> x;
 
  const FiniteElemSpace *feSpace = torus.getFESpace();
  const BasisFunction *basFcts = feSpace->getBasisFcts();
  int numBasFcts = basFcts->getNumber();
  DegreeOfFreedom *localIndices = GET_MEMORY(DegreeOfFreedom, numBasFcts);
  DOFAdmin *admin = feSpace->getAdmin();
 
  WorldVector<DOFVector<double>*> parametricCoords;
  for(i = 0; i < dow; i++) {
    parametricCoords[i] = NEW DOFVector<double>(feSpace, 
                                                "parametric coords");
  }
\end{lstlisting}

In the next step, we store the rotated vertex coordinates of the mesh in \verb+parametricCoords+,a vector of DOF vectors, where the first vector stores the first component of each vertex coordinate, and so on. In the STL map \verb+visited+, we store which vertices have already been visited, to avoid multiple rotations of the same point.

\begin{lstlisting}{}
  std::map<DegreeOfFreedom, bool> visited;
  TraverseStack stack;
  ElInfo *elInfo = stack.traverseFirst(torus.getMesh(), -1, 
				       Mesh::CALL_LEAF_EL | 
                                       Mesh::FILL_COORDS);
  while(elInfo) {
    basFcts->getLocalIndices(elInfo->getElement(), admin, localIndices);
    for(i = 0; i < dim + 1; i++) {
      dof = localIndices[i];
      x = elInfo->getCoord(i);
      YRotation::rotate(x, rotation);
      if(!visited[dof]) {
        for(j = 0; j < dow; j++) {
          (*(parametricCoords[j]))[dof] = x[j];
        }
        visited[dof] = true;
      }
    }
    elInfo = stack.traverseNext(elInfo);
  }
\end{lstlisting}

We create an instance of class \verb+ParametricFirstOrder+ which then is handed to the mesh. Now, in all future mesh traverses the vertex coordinates stored in \verb+parametricCoords+ are returned, instead of the original coordinates.

\begin{lstlisting}{}
  ParametricFirstOrder parametric(&parametricCoords);
  torus.getMesh()->setParametric(&parametric);
\end{lstlisting}

We rotate the right hand side function, reset \verb+adaptInfo+ and start the adaptation loop again. Now, we compute the solution on the rotated torus, which then is written to the files \verb+rotation1.mesh+ and \verb+rotation1.dat+.

\begin{lstlisting}{}
  f.rotate(rotation);
  adaptInfo->reset();
  adapt->adapt();

  DataCollector *dc = NEW DataCollector(feSpace, torus.getSolution());
  MacroWriter::writeMacro(dc, "output/rotation1.mesh");
  ValueWriter::writeValues(dc, "output/rotation1.dat");
  DELETE dc;
\end{lstlisting}

We perform another rotation. All we have to do is to modify the coordinates in \verb+parametricCoords+ and to inform $f$ about the rotation.

\begin{lstlisting}{}
  visited.clear();
  elInfo = stack.traverseFirst(torus.getMesh(), -1, 
			       Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
  while(elInfo) {
    basFcts->getLocalIndices(elInfo->getElement(), admin, localIndices);
    for(i = 0; i < dim + 1; i++) {
      dof = localIndices[i];
      x = elInfo->getCoord(i);
      YRotation::rotate(x, rotation);
      if(!visited[dof]) {
        for(j = 0; j < dow; j++) {
          (*(parametricCoords[j]))[dof] = x[j];
        }
        visited[dof] = true;
      }
    }
    elInfo = stack.traverseNext(elInfo);
  }

  f.rotate(rotation);
  adaptInfo->reset();
  adapt->adapt();

  dc = NEW DataCollector(feSpace, torus.getSolution());
  MacroWriter::writeMacro(dc, "output/rotation1.mesh");
  ValueWriter::writeValues(dc, "output/rotation1.dat");
  DELETE dc;
\end{lstlisting}

The solution is written to \verb+rotation2.mesh+ and \verb+rotation2.dat+.

Finally, we free some memory and finish the main program.

\begin{lstlisting}{}
  for(i = 0; i < dow; i++)
    DELETE parametricCoords[i];
  FREE_MEMORY(localIndices, DegreeOfFreedom, numBasFcts);
}
\end{lstlisting}

\subsection{Parameter file}

In the parameter file, we set the macro file to \verb+./macro/torus_macro.3d+. This two dimensional mesh is 8 times globally refined and successively projected on the torus.

\begin{lstlisting}{}
dimension of world:      3

torusMesh->macro file name:       ./macro/torus_macro.3d
torusMesh->global refinements:    8

torus->mesh:                 torusMesh
torus->dim:                  2
torus->polynomial degree:    1

torus->solver:                cg
torus->solver->max iteration: 1000
torus->solver->tolerance:     1.e-8
torus->solver->left precon:   diag
torus->estimator:             no
torus->marker:                no

torus->output->filename:       output/torus
torus->output->AMDiS format:   1
torus->output->AMDiS mesh ext: .mesh
torus->output->AMDiS data ext: .dat
\end{lstlisting}

\subsection{Macro file}

 The macro mesh defined in \verb+./macro/torus_macro.3d+ is shown in Figure \ref{f:torus macro}.
\begin{figure}
\center
\includegraphics[width=0.5\textwidth]{fig/torus_macro.jpg}
\caption{Macro mesh of the torus problem.}
\label{f:torus macro}
\end{figure}

\subsection{Output}

In Figure \ref{f:rotated torus}, the solutions on the three tori are shown.
\begin{figure}
\center
\includegraphics[width=0.7\textwidth]{fig/rotated_torus.jpg}
\caption{Solution on the original torus and on the two rotated tori.}
\label{f:rotated torus}
\end{figure}