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

leave a comment