[cig-commits] r11816 - cs/benchmark/cigma/trunk/src/tests
luis at geodynamics.org
luis at geodynamics.org
Tue Apr 15 08:10:04 PDT 2008
Author: luis
Date: 2008-04-15 08:10:03 -0700 (Tue, 15 Apr 2008)
New Revision: 11816
Added:
cs/benchmark/cigma/trunk/src/tests/TestLine.cpp
Log:
Compare function values along an arbitrary test line
Added: cs/benchmark/cigma/trunk/src/tests/TestLine.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/tests/TestLine.cpp (rev 0)
+++ cs/benchmark/cigma/trunk/src/tests/TestLine.cpp 2008-04-15 15:10:03 UTC (rev 11816)
@@ -0,0 +1,130 @@
+#include <iostream>
+#include <iomanip>
+#include <cmath>
+
+#include "../OkadaBenchmarkFields.h"
+#include "../FieldReader.h"
+#include "../PointField.h"
+
+using namespace std;
+using namespace cigma;
+
+
+double magnitude(double *y)
+{
+ return sqrt(y[0]*y[0] + y[1]*y[1] + y[2]*y[2]);
+}
+
+void print_point(double *y)
+{
+ cout << y[0] << " " << y[1] << " " << y[2] << " ";
+}
+
+int main()
+{
+ int N = 400;
+ double A[3] = { 0, 18000, -12000 };
+ double B[3] = { 24000, 18000, -12000 };
+
+ int i,j;
+ double t, dt;
+ dt = 1.0 / (N - 1);
+
+ double *x = new double[N*3];
+ for (i = 0; i < N; i++)
+ {
+ t = i * dt;
+ for (j = 0; j < 3; j++)
+ {
+ x[3*i + j] = A[j] + t * (B[j] - A[j]);
+ }
+ }
+
+ double *y[2];
+ y[0] = new double[N*3];
+ y[1] = new double[N*3];
+ for (i = 0; i < N; i++)
+ {
+ for (j = 0; j < 3; j++)
+ {
+ y[0][3*i + j] = 0.0;
+ y[1][3*i + j] = 0.0;
+ }
+ }
+
+
+ typedef okada::strike_slip::Disloc3d StrikeSlipDislocation;
+ StrikeSlipDislocation disloc3d;
+ for (i = 0; i < N; i++)
+ {
+ double *station = &x[3*i];
+ double *displacement = &y[0][3*i];
+ disloc3d.eval(station, displacement);
+ }
+
+
+ FieldReader fieldReader;
+ fieldReader.pointsReader.pointsPath = "../tet4_1000m.h5:/projection/points";
+ fieldReader.valuesReader.pointsPath = "../tet4_1000m.h5:/projection/data/snapshot0/displacements";
+ fieldReader.load_field();
+
+ Field *field = fieldReader.field;
+ if (field == 0)
+ {
+ cerr << "Could not load field!" << endl;
+ exit(1);
+ }
+ if (field->getType() != Field::POINT_FIELD)
+ {
+ cerr << "Expecting a point field" << endl;
+ exit(1);
+ }
+ //PointField *pointField = static_cast<PointField*>(fieldReader.field);
+ for (i = 0; i < N; i++)
+ {
+ double *station = &x[3*i];
+ double *displacement = &y[1][3*i];
+ field->eval(station, displacement);
+ }
+
+
+ double *residuals = new double[N*3];
+ for (i = 0; i < N; i++)
+ {
+ for (j = 0; j < 3; j++)
+ {
+ residuals[3*i + j] = y[1][3*i + j] - y[0][3*i + j];
+ }
+ }
+
+ cout << setprecision(12);
+ for (i = 0; i < N; i++)
+ {
+ t = i * dt;
+ cout << i << " ";
+ cout << t << " ";
+
+ double *X = &x[3*i];
+ double *Y0 = &y[0][3*i];
+ double *Y1 = &y[1][3*i];
+ double *D = &residuals[3*i];
+
+ print_point(X);
+ print_point(Y0);
+ print_point(Y1);
+ print_point(D);
+
+ cout << magnitude(Y0) << " ";
+ cout << magnitude(Y1) << " ";
+ cout << magnitude(D) << " ";
+ cout << endl;
+ }
+
+
+ delete [] x;
+ delete [] y[0];
+ delete [] y[1];
+ delete [] residuals;
+
+ return 0;
+}
More information about the cig-commits
mailing list