Mesh Mess

2D alpha shape demo in Mathematica

Posted in FEM, Numerical by Yi on 11/07/2009

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]]]

Tagged with: mathematica

leave a comment

Jacobi method demo

Posted in Numerical by Yi on 11/02/2009

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

Tagged with: algorithm, numerical

leave a comment

Seattle earthquake demo

Posted in Engineering by Yi on 10/27/2009

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.

Tagged with: demo, earthquake, video

leave a comment

ALE formulation

Posted in FEM by Yi on 10/26/2009

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 $\Omega$ and Lagrangian reference configuration domain $\Omega-0$ .

$x=\phi(\xi,t), \forall x\in \Omega,\xi\in \Omega_0$

Now assume a new reference region $\Omega_{\ast}$ , one have another map $f$ :

$x=f(c,t)=\phi(\xi,t),\forall c\in \Omega_{\ast}$

Here $f$ is identity mapping if $\Omega_{\ast}=\Omega$ while is $\phi$ if $\Omega_{\ast}=\Omega_0$ . Similar to material velocity in spatial configuration:

$u(x,t)=u(\phi(\xi,t),t)=\phi_{,t}(\xi,t)|_{\xi}$

new reference configuration also has a velocity, called “mesh velocity” in spatial configuration:

$u_{\ast}(x,t)=u_{\ast}(f(c,t),t)=f_{,t}(c,t)|_c$

With proper boundary condition $\phi(\xi,0)=\xi_0$ or $f(c,0)=c_0$ above two equations can be seen as PDE of $\phi$ or $f$ when corresponding velocities are given. Mesh motion $f$ is devised to reduce the mesh distortion during time evolution.

For a function $g(x),\forall x\in \Omega$ , there is counterpart of $g$ in other reference configurations:

$g(x,t)=g_0(\xi,t)=g_{\ast}(c,t)$

Take material derivatives:

$\frac{Dg(x,t)}{Dt}|_{\xi}=g_{,t}(x,t)+g_{,i}(x,t)x_{i,t}|_{\xi}=g_{\ast,t}(c,t)+g_{\ast,i}(c,t)c_{i,t}|_{\xi}$

in which $x_{i,t}|_{\xi}\equiv u_i$ is the velocity of material particle $\xi$ in spatial configuration, and $c_{i,t}|_{\xi}\equiv w_i$ is the velocity of material particle $\xi$ in arbitrary reference configuration. On the other hand, fixing material particle marker while taking time derivative of both sides of

$\phi_j(\xi,t)=f_j(c,t)$

gives

$u_j=f_{j,t}|_c+f_{j,i}c_{i,t}|_{\xi}=u_{\ast j}+x_{j,i}c_{i,t}|_{\xi}\Longrightarrow \frac{\partial x_j}{\partial c_i}w_i=u_j-u_{\ast j}$

where $u_{\ast}$ is the velocity of fixed mesh point $c$ in spatial configuration.

Combining last two conclusion we have

$Dg(x,t)/Dt=\dot{g_0}=g_{\ast,t}|_c+g_{,j}(u_j-u_{\ast j})$

Define $b\equiv u-u_{\ast}$ 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 $g$ and $g_{\ast}$ is

$Dg/Dt=g_{\ast,t}|_c+b\cdot \nabla_x g$

All the conservation laws in ALE formulation can be derived from expression above.

Tagged with: ALE, FEM, FSI

leave a comment

A sample of toy FDM code

Posted in Numerical by Yi on 10/22/2009

Solving equation

$u_t-u''=f(x,u)$

$f(x,u)=-2+2*(\sin(u)-\sin(x^2))$

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 $u(x)=x^2$ . Here is time evolution plot of the numerical solution.

Tagged with: Finite difference, numerical, time integration

1 comment

Older Posts »

M	T	W	T	F	S	S
« Oct
						1
2	3	4	5	6	7	8
9	10	11	12	13	14	15
16	17	18	19	20	21	22
23	24	25	26	27	28	29
30

Mesh Mess

2D alpha shape demo in Mathematica

Jacobi method demo

Seattle earthquake demo

ALE formulation

A sample of toy FDM code

Recent Posts

Tags

Meta