FEM, CFD, CAE, and all that.

## Tag: algorithm

### POD snippet in R

I’ve been working on proper orthogonal decomposition (POD) for a while. Though it’s nice to have packages such as modred out there, I cannot help myself writing a snippet using R, since lots of my previosu data processing was in that language. Without all the theory chitchat  like K-L transformation, from the the implementation-on-discrete-data perspective the (snapshot) POD method boils down to just singular value decomposition (SVD), and it’s out of box in R, which makes the task mostly straightforward. That’s why a POD function in R essentially calls svd in R (and it further calls LAPACK routines DGESDD and ZGESVD.):

pod <- function(arr, nmod=dim(arr)[1], nv=nmod){
pod <- svd(arr, nmod, nv)
return(pod)
}

Here “nmod” is the number of POD modes to be retrieved. Though it’s default to be the length of each snapshot, it’s rarely necessary this large. Mostly of the time a reasonable correlated set of snapshots requires less than 10 modes to obtain most of the system energy.

With pod of data “arr” obtained, we can do several thing immediately. We  can obtain the modes by retrieving the right singular vector:

pod_u<- function(pod, nnod){
u <- matrix(unlist(pod["u"]), nrow = nnod )
return(u)
}

or get the singular values:

pod_sv<- function(pod){
sv <- as.vector(unlist(pod["d"]))
return(sv)
}

or the coefficients of the pod modes for the approximation:

pod_coef<- function(pod, nshot, nmod){
v <- matrix(unlist(pod["v"]), nrow = nshot )
sv <- pod_sv(pod)[1:nmod]
coef <- diag(sv) %*% t(v)
return(coef)
}

and eventually the approximation itself

pod_approx<- function(pod, nshot, nmod, nnod){
u <- pod_u(pod, nnod)
coef <- pod_coef(pod, nshot, nmod)
approx <- u %*% coef
return(approx)
}

### Jacobi method demo

Code below is to solve the static version of the equation in previous FDM code demo, with B.C. changed as homogeneous. The Jacobi method is used to replace the LAPACK solver which gives inexact method. The exact solution is
$u(x)=-x^2/2+x/2$
Using initial guess as
$u_0(x)=\sin(2\pi x)+0.5\sin(10\pi x)+0.3\sin(20\pi x)+0.1\sin(40\pi x)$
one can see that most of the computing cost is spent on smoothing lowest mode.


#include <stdio.h>
#include <math.h>
#include <iostream>
#include <fstream>
//#include <vector>
//#include <sunperf.h>
using namespace std;
//typedef std::vector<float> fvector_s;

//Maximum of an array.
float max(float a[], int n) {
float v = a[0];
for (int i=1; i<n; i++) {
if (a[i] > v) {
v = a[i];
}
}
return v;
}

int main ()
{
int n=500;
float tmp, b[n], a_d[n], a_u[n-1], a_l[n-1],    \
tol=1.e-7, x0[n], x1[n], x2[n], error=1.;
float h=1./(n+1), pi=3.1416;
int count=0;

//initializing a
for(int i=0; i<n; ++i){
a_d[i]=2.0;
}
for(int i=0; i<n-1; i++){
a_l[i]=-1.;
a_u[i]=-1.;
}

//initializaing b[]
for(int i=0; i<n; ++i){
b[i]=h*h;
}

//initial guess
for(int i=0; i<n; ++i){
x0[i]=sin(2*pi*h*(i+1))+0.5*sin(10*pi*h*(i+1))+\
0.3*sin(20*pi*h*(i+1))+0.1*sin(40*pi*h*(i+1));
x1[i]=x0[i];
}

ofstream outfile;
outfile.open("jac.out");
do{
count+=1;
for(int i=0; i<n; ++i){
if(i==0)
x2[i]=(b[i]-a_u[i]*x1[i+1])/a_d[i];
else if(i==n-1)
x2[i]=(b[i]-a_l[i-1]*x1[i-1])/a_d[i];
else
x2[i]=(b[i]-a_u[i]*x1[i+1]-a_l[i-1]*x1[i-1])/a_d[i];
}
for(int i=0; i<n; i++){
x1[i]=x2[i]-x1[i];
}
error=max(x1, n)/max(x0, n);
for(int i=0; i<n; i++){
x1[i]=fabs(x2[i]-(-(i*h)*(i*h)/2+i*h/2));
}
outfile << count << "\t" << error << endl;
for(int i=0; i<n; ++i){
x1[i]=x2[i];
}
}while(error>tol);
cout << "Number of iterations:" << count << endl;
for(int i=0; i<n; i++){
outfile << count << "\t" << x2[i] << endl;
}

outfile.close();
return 0;
}


Result of Jacobi method at different iteration

### Alpha shape

Alpha shape is a generalized way determining the geometric objects by a set of points, of which convex hull is a special case. It can be applied to FEM for boundary recognition. Essentially this method identifies the boundary of a set of points by eliminating those triangles (generated by triangulation of the points) which can be put into empty circles whose diagram are defined by threshold alpha. It is intuitive for a more simple criteria: eliminating the triangles whose circumcircles contains no points in the defined set. Below is a simple demo by Mathematica which tries to identify the boundary of a subset of points within circle radius 5 from a set points randomly distributed. The red arc isolates the point set $P$  inside, the first step is to eliminate the triangles whose circumcirlces contains no points in $P$, the next step is to eliminate the triangles contain nodes other than those in $P$. More delicate results can be reached by tuning alpha circle instead of using simple circumcirlces.

### Bowyer/Watson algorithm: 3D

3D version implementation of Bowyer/Watson algorithm is almost the same as the 2D. Instead of eliminating lines that belong only to the triangles whose circumcirlce contains the newly added point, in 3D we eliminate the facets that belong only to the tetrahedra whose circumsphere contains the newly added point, then form new tetrahedra by appending new point to those facets’ nodes. Here is an example generated by Mathematica scripts I compiled. Below on the right is tetrahedronization of 200 random points, together with enclosing tetrahedron. On the left is he tetrahedrons generated by first two points.

### DistMesh

DistMesh is no doubt the most interesting way meshing a domain, simple and elegant. There is a MATHEMATICA version of the original Matlab code, but only in 2D. Base on that, I wrote a script for mesh test. Hopefully continue working on this one would bring some application to my research. Here is one simple case of generating a non-uniform mesh within a circle, whose distance funtion is the easist to define. At the end of termination, most of the ill-shaped triangles are corrected.