# Discrete K-surfaces tutorial

This tutorial explains how one may implement the discrete K-surfaces application with jReality. This tutorial is one of the longer tutorials which explains a small but still reasonably complete application. The formulas behind this application are explained in Ulrich Pinkall: Designing Cylinders with Constant Negative Curvature.

## Vector operations in $R^3$

For the vector calculus in $R^3$ we represent vectors as arrays of doubles of length 3 and implement operations on these vectors as static methods of a class named R3:


/** This class provides simple static methods to do verctor calculation in R<sup>3</sup>.
* The vectors of R<sup>3</sup> are represented by array of 3 doubles.
*
*/
public class R3 {
/** target = a + b
*
* @param a array of 3 doubles
* @param b array of 3 doubles
* @param target array of 3 doubles
*/
public final static void plus(double[] a, double[] b, double[] target) {
target[0] = a[0] + b[0];
target[1] = a[1] + b[1];
target[2] = a[2] + b[2];
}

/** target = a-b
*
* @param a array of 3 doubles
* @param b array of 3 doubles
* @param target array of 3 doubles
*/
public final static void minus(double[] a, double[] b, double[] target) {
target[0] = a[0] - b[0];
target[1] = a[1] - b[1];
target[2] = a[2] - b[2];
}

/** target = a x b, i.e., the <a href="http://en.wikipedia.org/wiki/Crossproduct">cross product or vector product</a><br>
*
* WARNING: This method fails when the target is same as one of its arguments.
*
* @param a array of 3 doubles
* @param b array of 3 doubles
* @param target array of 3 doubles
*/
public final static void cross(double[] a, double[] b, double[] target) {
double x = -a[1]*b[2] + a[2]*b[1];
double y = -a[2]*b[0] + a[0]*b[2];
target[2] = -a[0]*b[1] + a[1]*b[0];
target[0] = x;
target[1] = y;
}

/** target = a, copy the entries of a into target
*
* @param a array of 3 doubles
* @param target array of 3 doubles
*/
public final static void copy(double[] a, double[] target) {
target[0] = a[0];
target[1] = a[1];
target[2] = a[2];
}

/** target = 0, reset the three entries of target to 0
*
* @param target array of 3 doubles
*/
public final static void zero(double[] target) {
target[0] = 0;
target[1] = 0;
target[2] = 0;
}

/** &lt;a,b&gt;, i.e., the <a href="http://en.wikipedia.org/wiki/Dotproduct">dot product or inner product</a>
*
* @param a array of 3 doubles
* @param b array of 3 doubles
* @return the value of the dot product
*/
public final static double dot(double[] a, double[] b) {
return a[0]*b[0] + a[1]*b[1] + a[2]*b[2];
}

/** target = s * a, product of a real with a vector
*
* @param s the scalar real
* @param a array of 3 doubles
* @param target array of 3 doubles
*/
public final static void times(double s, double[]a, double[]target) {
target[0] = s * a[0];
target[1] = s * a[1];
target[2] = s * a[2];
}

/** |a|, the length of a vector.
*
* @param a array of 3 doubles
* @return the length of a
*/
public static double norm(double[] a) {
return Math.sqrt(normSquared(a));
}

/** |a|^2, the square of the length of a vector.
*
* @param a array of 3 doubles
* @return the square of the length of a
*/
public static double normSquared(double[] a) {
return a[0]*a[0]+a[1]*a[1]+a[2]*a[2];
}
}

## Calculate the Gauss map

The Gaus map of a discrete K-surface consists of spherical parallelograms, i.e., spherical quadrilateral that are invariant under a 180° rotation. But the parallelograms are skew in our coordinates (this corresponds to the coordinate change $u=x+t,\ v=x-t$ which transforms the form $N_{xx}-N_{tt}=0$ of the wave equation into $N_{uv}=0$. The equation $f_{uv}=0$), i.e., a parallelogram is formed by the Gauss map at (i-2,j-1), (i-1,j-1), (i,j), and (i-1,j). Thus the Gauss map at (i,j) is 180° rotation of the Gauss map at (i-2,j-1) about the mid point of the diagonal between (i-1,j-1) and (i-1,j).

Let us implement the algorithm for computing the Gauss map as a static method of the class KSurfaces:

public class KSurfaces {
/** Computes the Gauss map of a K-surface from initial Cauchy data, i.e., two closed curves. Assums that the points of
* the initial data lie on the 2-sphere this method generates more points so that the quadrilaterals
* (i-1,j-1), (i-1,j), (i,j),(i-2,j-1) form spherical parallelograms.
*
* @param initialAnnulus the initial data, i.e., a double array double[2][n][3] containing 2 polygons with n vertices,
*  which are both interpreted as closed curves connecting that last and the first one.
* @param target a double array that may be filled with the result, i.e., an array double[m][n][3], where m&gt;1 is the number of
* 	time steps to be calculated.
*/
public static void gaussMapFromInitialAnnulus(double[][][] initialAnnulus, double[][][] target) {
final int m = target.length;
final int n = target[0].length;
final double[] a = new double[3];

// copy first two rows
for (int i=0; i<2; i++) {
for (int j=0; j<n; j++) {
R3.copy(initialAnnulus[i][j], target[i][j]);
}
}

// compute the other rows
for (int i=2; i<m; i++) {
for (int j=0; j<n; j++) {
int k = j==0 ? n-1 : j-1;
R3.plus(target[i-1][k], target[i-1][j], a);
double[] w = target[i-2][k];
double s = 2 * R3.dot(a,w) / R3.dot(a,a);
R3.times(s, a, a);
R3.minus(a, w, target[i][j]);
}
}
}
}

## Display the Gauss map

Now we are able to display the Gauss map. For the moment, we choose an ad hoc initial condition: the intersection of two elliptic cones around the z-achses with the 2-sphere. Lets put the code for this in yet another class named KSurfacesApp.

### The fields of KSurfacesApp

import de.jreality.geometry.IndexedFaceSetFactory;
import de.jreality.plugin.JRViewer;
import de.jreality.plugin.JRViewer.ContentType;
import de.jreality.plugin.basic.ViewShrinkPanelPlugin;
import de.jreality.plugin.content.ContentAppearance;
import de.jreality.plugin.content.ContentTools;
import de.jreality.scene.SceneGraphComponent;
import de.jtem.beans.InspectorPanel;
import de.varylab.jrworkspace.plugin.PluginInfo;

/** This class contains the application to display the K-surfaces and their Gauss map.
*/
public class KSurfacesApp {
private int m=3;
private int n=30;
private double a1=1;
private double a2=.9;
private double b1=1.2;
private double b2=1.1;
private double[][][] gaussMap;
private double[][][] initialGaussMap;

final private SceneGraphComponent root=new SceneGraphComponent("k-surfaces");
final private SceneGraphComponent gaussMapSgc=new SceneGraphComponent("Gauss map");
final private IndexedFaceSetFactory gaussMapIfs = new IndexedFaceSetFactory();

### The constructor of KSurfacesApp

public KSurfacesApp() {
gaussMapIfs.setGenerateEdgesFromFaces(true);
gaussMapIfs.setGenerateFaceNormals(true);
gaussMapSgc.setGeometry(gaussMapIfs.getGeometry());

update();
}

### The method that generates the initial condition

This method generates two generations of normal vectors for use in the method gaussMapFromInitialAnnulus

private void ellipticConesInitialCondition() {
initialGaussMap = new double[2][n][3];
for (int j=0; j<n; j++){
double s=a1*Math.sin(2*Math.PI/n*j);
double c=a2*Math.cos(2*Math.PI/n*j);
double f=1/Math.sqrt(s*s+c*c+1);
initialGaussMap[0][j][0]=s*f;
initialGaussMap[0][j][1]=c*f;
initialGaussMap[0][j][2]=1*f;

s=b1*Math.sin(2*Math.PI/n*(j-.5));
c=b2*Math.cos(2*Math.PI/n*(j-.5));
f=1/Math.sqrt(s*s+c*c+1);
initialGaussMap[1][j][0]=s*f;
initialGaussMap[1][j][1]=c*f;
initialGaussMap[1][j][2]=1*f;
}
}

### Recalculate the gauss map and update the display

private void gaussMap() {
//calculate the gauss map
gaussMap=new double[m][n][3];
KSurfaces.gaussMapFromInitialAnnulus(initialGaussMap, gaussMap);

//rearrange the vertices for use in an IndexedFaceSetFactory
double[][] vertices=new double[n*m][3];
for (int i = 0; i < m; i++) {
for (int j = 0; j < n; j++) {
for (int k = 0; k < 3; k++) {
vertices[j+i*n][k]=gaussMap[i][j][k];
}
}
}

//calculate the combinatorics of the faces for use in an IndexedFaceSetFactory
int[][] faces = new int[n*(m-2)][4];
for (int i = 0; i < m-2; i++) {
for (int j = 0; j < n; j++) {
int k = j + i*n;
faces[k][0] = j       + i*n;
faces[k][1] = j       + (i+1)*n;
faces[k][2] = (j+1)%n + (i+2)*n;
faces[k][3] = (j+1)%n + (i+1)*n;
}
}
//put the data into the IndexedFaceSetFactory ifs and update the Geometry
gaussMapIfs.setVertexCount(n*m);
gaussMapIfs.setVertexCoordinates(vertices);
gaussMapIfs.setFaceCount(n*(m-2));
gaussMapIfs.setFaceIndices(faces);
gaussMapIfs.update();
}

### Getters and Setters


public int getM() { return m; }
public void setM(int m) {
if (m<2) return;
this.m = m;
}

public int getN() { return n; }
public void setN(int n) {
if (n<2) return;
this.n = n;
}

public double getA1() {	return a1; }
public void setA1(double a1) {	this.a1 = a1; 	}

public double getA2() { return a2; }
public void setA2(double a2) { this.a2 = a2; }

public double getB1() {	return b1; }
public void setB1(double b1) {	this.b1 = b1; }

public double getB2() {	return b2; }
public void setB2(double b2) { this.b2 = b2; }

/**
* @return the root of the scene graph with the Gauss map of the k-surface
*/
public SceneGraphComponent getRoot() {
return root;
}

### Update method

/** This method needs to be called whenever the changes of the parameters should have an effect on the
* scene graph one gets from {@link #getRoot()}.
*/
public void update() {
ellipticConesInitialCondition();
gaussMap();
}

### The main method

public static void main(String[] args) {
KSurfacesApp ksa = new KSurfacesApp();

//initialize the ViewerVR
JRViewer v = new JRViewer();
v.registerPlugin(new ContentAppearance());
v.registerPlugin(new ContentTools());
v.setContent(ksa.getRoot());

//add an Inspector for the KSurfacesApp to the viewer
InspectorPanel inspector = new InspectorPanel();
inspector.setObject(ksa, "update");
ViewShrinkPanelPlugin plugin = new ViewShrinkPanelPlugin() {
public PluginInfo getPluginInfo() {
return new PluginInfo("Domain");
}
};
v.registerPlugin(plugin);

v.startup();
}

## Editor for the initial curve and velocity

Have a look at the Polygon editor tutorial and use this for the KSurfacesApp class.

Here is the essential code you may add to the constructor

//initial condition at time 1 draggable points
a1=0;a2=0;
b1=.2;b2=.2;
ellipticConesInitialCondition();
//to keep the points on the unit sphre override the transform method
dragPointSet=new DragPointSet(initialGaussMap[1]) {
@Override
public void transform(double[] vertex) {
double l=R3.norm(vertex);
vertex[0]=vertex[0]/l;
vertex[1]=vertex[1]/l;
vertex[2]=vertex[2]/l;
}
};
public void stateChanged(ChangeEvent e) {
KSurfacesApp.this.update();
}
});

//initial condition at time 1 subdivided Polygon points
time1SubdivPoly = new SubdividedPolygon(dragPointSet);
time1SubdivPoly.setSubdivisionLevel(0);
PointSequenceView time1CurveView = new PointSequenceView(time1SubdivPoly);
time1CurveView.setPointColor(Color.yellow);
time1CurveView.setLineColor(Color.green);
gaussMapSgc.addChild(time1CurveView.getBase());

And two lines two the update() method

n=time1SubdivPoly.getPoints().length;
ellipticConesInitialCondition();
initialGaussMap[1]=time1SubdivPoly.getPoints();

## Integrate the Gauss map to get the surface

The first row, i.e., to calculate f(0,j) we use the formula $f_{\tilde d}=f_d +(N_d+N_{\tilde d})\times N_{r}$, where in our xt-coordinates $d=(0,j),\ \tilde d=(0,j+1),\ r=(1,j+1)$. For the other rows we use $f_r=f_d+N_d\times N_r$, where $d=(i-1,j),\ r=(i,j)$. Let us put this algorithm as a static method into the class KSurfaces:

	/** Calculates the K-surface from the given Gauss map, assuming that the Gauss map consists
* of spherical parallelograms as described in {@link #gaussMapFromInitialAnnulus(double[][][], double[][][])}.
*
* @param gaussMap the given Gauss map.
* @param target a double array with enough space to hold the resulting surface. A quad mesh where the
* quadrilaterals consist of (i-1,j-1), (i-1,j), (i,j),(i-2,j-1).
*/
public static void kSurfaceFromGaussMap(double[][][] gaussMap, double[][][] target) {
final int m = gaussMap.length;
final int n = gaussMap[0].length;
final double[] a = new double[3];

// compute first row
R3.zero(target[0][0]);
for (int j=0; j<n-1; j++) {
int k = (j+1)%n;
R3.plus(gaussMap[0][j], gaussMap[0][k], a);
R3.cross(gaussMap[1][k], a, a);
R3.plus(target[0][j], a, target[0][k]);
}

// compute the other rows
for (int i=1; i<m; i++) {
for (int j=0; j<n; j++) {
R3.cross(gaussMap[i-1][j], gaussMap[i][j], a);
R3.plus(target[i-1][j], a, target[i][j]);
}
}
}

## Display the surface

We add to the class KSurfacesApp two new imports

import de.jreality.math.MatrixBuilder;
import de.jreality.scene.Transformation; 

Then a field for the map, its scene graph component and its factory:

private double[][][] map;
final private SceneGraphComponent mapSgc=new SceneGraphComponent("The K-surface");
final private IndexedFaceSetFactory mapIfs = new IndexedFaceSetFactory();

Then the following code to the constructor KSurfacesApp()

mapIfs.setGenerateEdgesFromFaces(true);
mapIfs.setGenerateFaceNormals(true);
mapSgc.setGeometry(mapIfs.getGeometry());
mapSgc.setTransformation(new Transformation(
MatrixBuilder.euclidean().translate(3,0,0).getArray()));
root.addChild(mapSgc);

The we refactor the methode gaussMap() to avoid code doubling. Just replace the old version by

//calculate the gauss map and update it via gaussMapIfs
private void gaussMap() {
//calculate the gauss map
gaussMap=new double[m][n][3];
KSurfaces.gaussMapFromInitialAnnulus(initialGaussMap, gaussMap);

double[][] vertices = rearrangeVertices(gaussMap);
int[][] faces = calcFaces();

//put the data into the IndexedFaceSetFactory and update the Geometry
gaussMapIfs.setVertexCount(n*m);
gaussMapIfs.setVertexCoordinates(vertices);
gaussMapIfs.setFaceCount(n*(m-2));
gaussMapIfs.setFaceIndices(faces);
gaussMapIfs.update();
}

//calculate the map of K-surface and update it via gaussMapIfs
private void map() {
//calculate the map
map=new double[m][n][3];
KSurfaces.kSurfaceFromGaussMap(gaussMap, map);

double[][] vertices = rearrangeVertices(map);
int[][] faces = calcFaces();

//put the data into the IndexedFaceSetFactory and update the Geometry
mapIfs.setVertexCount(n*m);
mapIfs.setVertexCoordinates(vertices);
mapIfs.setFaceCount(n*(m-2));
mapIfs.setFaceIndices(faces);
mapIfs.update();
}

//calculate the combinatorics of the faces for use in an IndexedFaceSetFactory
private int[][] calcFaces() {
int[][] faces = new int[n*(m-2)][4];
for (int i = 0; i < m-2; i++) {
for (int j = 0; j < n; j++) {
int k = j + i*n;
faces[k][0] = j       + i*n;
faces[k][1] = j       + (i+1)*n;
faces[k][2] = (j+1)%n + (i+2)*n;
faces[k][3] = (j+1)%n + (i+1)*n;
}
}
return faces;
}

//rearrange the vertices for use in an IndexedFaceSetFactory
private double[][] rearrangeVertices(double[][][] xtVertices) {
double[][] vertices=new double[n*m][3];
for (int i = 0; i < m; i++) {
for (int j = 0; j < n; j++) {
for (int k = 0; k < 3; k++) {
vertices[j+i*n][k]=xtVertices[i][j][k];
}
}
}
return vertices;
}

Finally call the method map() from update()

map();