[cig-commits] r13860 - cs/cigma/trunk/src

[luis at geodynamics.org]

Author: luis
Date: 2009-01-13 03:12:04 -0800 (Tue, 13 Jan 2009)
New Revision: 13860

Added:
   cs/cigma/trunk/src/fn_disloc3d.cpp
   cs/cigma/trunk/src/fn_disloc3d.h
Log:
OkadaDisloc3d function that we can register later on.

This is a nontrivial example that demonstrates how to write
your own analytic function by calling an "external" wrapper.

Added: cs/cigma/trunk/src/fn_disloc3d.cpp
===================================================================
--- cs/cigma/trunk/src/fn_disloc3d.cpp	                        (rev 0)
+++ cs/cigma/trunk/src/fn_disloc3d.cpp	2009-01-13 11:12:04 UTC (rev 13860)
@@ -0,0 +1,101 @@
+#include "fn_disloc3d.h"
+#include "core_misc.h"
+#include "cervelli.h"
+#include "disloc3d.h"
+
+using namespace std;
+using namespace benchmark::quasistatic;
+
+strikeslipnog::OkadaDisloc3d::OkadaDisloc3d()
+{
+    // elastic parameters
+    mu = 0.25;
+    nu = 0.25;
+
+    // allocate fault info
+    nsubfaults = 400;
+    subfaults.resize(nsubfaults * 11);
+    models.resize(nsubfaults * 10);
+    flag2.resize(nsubfaults);
+
+    // define strike-slip fault
+    std::valarray<double> n;
+
+    double taperd = (16000.0 - 12000.0)/400;  // 10.0
+    linspace(n, 12000.0, 16000.0, nsubfaults);
+
+    double *fi;
+    double fx1, fy1;
+    double fx2, fy2;
+    double dip, D, Pr, bd;
+    double ss, ds, ts;
+
+    for (i = 0; i < nsubfaults; i++)
+    {
+        fi = &subfaults[11*i];
+
+        fx1 = 12000.0;
+        fy1 = -n[i];
+
+        fx2 = 12000.0;
+        fy2 = +n[i];
+
+        dip = 90.0;
+        D   = n[i];
+        Pr  = 0.25;
+        bd  = 0;
+        ss  = -1.0/nsubfaults;
+        ds  = 0;
+        ts  = 0;
+
+        fi[0]  = fx1;
+        fi[1]  = fy1;
+        fi[2]  = fx2;
+        fi[3]  = fy2;
+        fi[4]  = dip;
+        fi[5]  = D;
+        fi[6]  = Pr;
+        fi[7]  = bd;
+        fi[8]  = ss;
+        fi[9]  = ds;
+        fi[10] = ts;
+    }
+
+    // Get the "model matrix" with fault patch info
+    // in the form for cervelli calcs
+    getM(&subfaults[0], nsubfaults, mu, &models[0]);
+}
+
+strikeslipnog::OkadaDisloc3d::~OkadaDisloc3d()
+{
+}
+
+bool strikeslipnog::OkadaDisloc3d::eval(double *station, double *disloc)
+{
+    double flag = -1;
+
+    // Call the wrapper function
+    disloc3d(&models[0], nsubfaults,
+             station, 1,
+             mu, nu,
+             displacement, derivatives, stress,
+             &flag, &flag2[0]);
+
+    if (flag > 0)
+    {
+        return false;
+    }
+
+    // copy answer to the output array
+    disloc[0] = displacement[0];
+    disloc[1] = displacement[1];
+    disloc[2] = displacement[2];
+
+    //
+    // Note: how should we select other fields? access directly perhaps.
+    // Value should be set correctly after call to disloc3d wrapper.
+    //
+
+    return true;
+}
+

Added: cs/cigma/trunk/src/fn_disloc3d.h
===================================================================
--- cs/cigma/trunk/src/fn_disloc3d.h	                        (rev 0)
+++ cs/cigma/trunk/src/fn_disloc3d.h	2009-01-13 11:12:04 UTC (rev 13860)
@@ -0,0 +1,49 @@
+#ifndef FN_DISLOC3D_H
+#define FN_DISLOC3D_H
+
+#include "Function.h"
+#include <valarray>
+
+namespace benchmark
+{
+    namespace quasistatic
+    {
+        namespace strikeslipnog
+        {
+            class OkadaDisloc3d;
+        }
+
+        //namespace reversenog
+        //{
+        //    class OkadaDisloc3d;
+        //}
+    }
+}
+
+
+class benchmark::quasistatic::strikeslipnog::OkadaDisloc3d : public cigma::Function
+{
+public:
+    OkadaDisloc3d();
+    ~OkadaDisloc3d();
+
+    virtual int n_dim() { return 3; }
+    virtual int n_rank() { return 3; }
+
+    virtual bool eval(double *station, double *disloc);
+
+public:
+    double mu, nu;
+    double displacement[3];
+    double derivatives[9];
+    double stress[6];
+
+public:
+    int nsubfaults;
+    std::valarray<double> subfaults;
+    std::valarray<double> models;
+    std::valarray<double> flag2;
+};
+
+
+#endif