2D alpha shape demo in Mathematica
This is the source code for a past post.
<< ComputationalGeometry`;
(*Loading test points*)
p = RandomReal[{0, 10}, {1000, 2}];
p = Select[p, Norm[# - {5, 5}] < 5 &]; p = Select[p, Norm[# - {7.5, 5}] > 2.5 &];
ListPlot[p, PlotRange -> Full, AspectRatio -> 1]
tri = DelaunayTriangulation[p];
tri = Thread[List @@ #] & /@ tri // Flatten[#, 1] &;
tri = Union[Sort /@ tri];
(*Main function def*)
(*0 for interior bar,1 for boundary bar, 2 for bar to be erased*)
ClearAll[f];
emptyQ[center_, alpha_, plist_] := Module[
{empty = True, n = 1},
While[empty && n <= Length[plist], empty = Norm[plist[[n]] - center] > alpha; n++];
empty
];
f[alpha_, plist_, {id1_, id2_}] := Module[
{p1 = plist[[id1]], p2 = plist[[id2]], center1, center2, lhalf,
center1emptyQ, center2emptyQ, property},
lhalf = Norm[p2 - p1]/2;
center1 = (p2 + p1)/2 +
Sqrt[(alpha/lhalf)^2 - 1] {{0, -1}, {1, 0}}.((p2 - p1)/2);
center2 = (p2 + p1)/2 +
Sqrt[(alpha/lhalf)^2 - 1] {{0, 1}, {-1, 0}}.((p2 - p1)/2);
center1emptyQ =
emptyQ[center1, alpha, Delete[plist, {{id1}, {id2}}]];
center2emptyQ =
emptyQ[center2, alpha, Delete[plist, {{id1}, {id2}}]];
Which[center1emptyQ \[And] center2emptyQ,
2, (! center1emptyQ) \[And] (! center2emptyQ),
0, (center1emptyQ \[And] (! center2emptyQ)) \[Or] ((!
center1emptyQ) \[And] center2emptyQ), 1]
(*Print[property];
Show[Graphics[{Red,Circle[center1,alpha],Circle[center2,alpha]}],
ListPlot[plist,AspectRatio->1]]*)
] /; alpha > Norm[plist[[id1]] - plist[[id2]]]/2
f[alpha_, plist_, {id1_, id2_}] :=
2 /; alpha < Norm[plist[[id1]] - plist[[id2]]]/2
(*plot a test result*)
test = Extract[tri, Position[f[0.33, p, #] & /@ tri, 1]];
Graphics[GraphicsComplex[p, Line[test]]]
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
Using initial guess as
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;
}
Seattle earthquake demo
From today’s post of Dave’s landslide blog, this is quite some nice demo of seawall failure and viaduct collapse assuming an long duration earthquake attack at Seattle from sea.
ALE formulation
There is a old post for Lagrangian configuration, which is popular in solid mechanics but not in CFD. ALE is dedicated for combining both Eularian and Lagrangian configuration. In Lagrangian, one has the map between Eularian configuration region and Lagrangian reference configuration domain
.
Now assume a new reference region , one have another map
:
Here is identity mapping if
while is
if
. Similar to material velocity in spatial configuration:
new reference configuration also has a velocity, called “mesh velocity” in spatial configuration:
With proper boundary condition or
above two equations can be seen as PDE of
or
when corresponding velocities are given. Mesh motion
is devised to reduce the mesh distortion during time evolution.
For a function , there is counterpart of
in other reference configurations:
Take material derivatives:
in which is the velocity of material particle
in spatial configuration, and
is the velocity of material particle
in arbitrary reference configuration. On the other hand, fixing material particle marker while taking time derivative of both sides of
gives
where is the velocity of fixed mesh point
in spatial configuration.
Combining last two conclusion we have
Define as the difference between material velocity and mesh velocity in spatial configuration, when both material particle marker and mesh point marker are at the same spatial position at that moment. Finally, the transformation of
and
is
All the conservation laws in ALE formulation can be derived from expression above.
A sample of toy FDM code
Solving equation
Time integration is full implicit, and during each time step Newton’s method is employed.
PROGRAM MAIN
!--------------------------------------------------
!Implicit time evolution of equation
!u_t-u''=f(x,u)
!f(x,u)=-2+r*(sin(u)-sin(x^2)),r=2
!--------------------------------------------------
! t=time points
! dt=time step
! h=neighboring nodes distance
! u=unknown
! x=node location
! n=No. of unknowns
!USE SUNPERF
IMPLICIT NONE
INTEGER(4) :: n
REAL(4) :: elapse,time(2),etime
REAL(8) :: t=0,dt=0.1
REAL(8) :: h
REAL(8),ALLOCATABLE :: u(:),x(:)
INTEGER(4) :: count,i
n=50
h=1/(DBLE(n)+1.)
ALLOCATE (x(n),u(n))
x=(/(i*h,i=1,n)/)
u=(/((i*h)**0.5,i=1,n)/)
OPEN(UNIT=100,FILE='data',STATUS='UNKNOWN',ACTION='&
&WRITE')
DO i=1,n
WRITE(100,*) t,x(i),u(i)
END DO
DO WHILE (t<1.0)
CALL newton(n,h,dt,u,x)
t=t+dt
DO i=1,n
WRITE(100,*) t,x(i),u(i)
END DO END DO CLOSE(UNIT=100)
elapse=etime(time)
WRITE(*,*) 'elapsed:',elapse
END PROGRAM MAIN
!--------------------------------------------------
SUBROUTINE newton(n,h,dt,u,x)
IMPLICIT NONE
INTEGER(4),INTENT(IN) :: n
REAL(8),INTENT(IN) :: x(n),h,dt
REAL(8),INTENT(INOUT) :: u(n)
REAL(8) :: r=2,c(n),error,b(n,1),b0(n,1),a(2,n)
REAL(8) :: alpha=1.0,beta=-1.0
INTEGER(4) :: info,i,iter
c=u*h*h/dt
iter=0
error=1
DO WHILE (error>1.e-13)
iter=iter+1
a(1,:)=-1
a(2,:)=2+h*h/dt
b(:,1)=(-2+r*(SIN(u)-SIN(x*x)))*h*h+c
b(n,1)=b(n,1)+1
CALL DSBMV('u',n,1,alpha,a,2,u,1,beta,b,1)
b=-b
a(2,:)=a(2,:)-r*COS(u)*h*h
CALL DGTSV(n,1,a(1,2:),a(2,:),a(1,2:),b,n,info)
error=MAXVAL(ABS(b))
u=u+b(:,1)
END DO
WRITE(*,*) 'INFO',info,'Number of iterations:',iter
END SUBROUTINE newton
!--------------------------------------------------
Output would be like:
% f95 -dalign implicit_newton.f95 -xlic_lib=sunperf % a.out INFO 0 Number of iterations: 4 INFO 0 Number of iterations: 4 INFO 0 Number of iterations: 4 INFO 0 Number of iterations: 4 INFO 0 Number of iterations: 4 INFO 0 Number of iterations: 3 INFO 0 Number of iterations: 3 INFO 0 Number of iterations: 3 INFO 0 Number of iterations: 3 INFO 0 Number of iterations: 3 elapsed: 0.061407
Of course the exact solution is . Here is time evolution plot of the numerical solution.

