Secrete of a craftsman

This post in SA about Antonio Stradivari’s secrete of violin building makes great sense in other crafts work as well: it is not the tool that defines a master, it is the skill. For my own experience,  once for a long time I have been obsessed in looking for a perfect tool in everything, especially in work on coding, on computing, among others. I was told finally, in particular in scientific computing business,  you do not choose your tools, your tools choose you.

Istanbul Earthquake-Safe Building

This is awesome. I suppose the interactive behavior of all those foundations is a large issue.

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;
}

Result of Jacobi method at different iteration

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.