[cig-commits] r11856 - cs/benchmark/cigma/trunk/src/tests
luis at geodynamics.org
luis at geodynamics.org
Wed Apr 23 09:56:43 PDT 2008
Author: luis
Date: 2008-04-23 09:56:43 -0700 (Wed, 23 Apr 2008)
New Revision: 11856
Added:
cs/benchmark/cigma/trunk/src/tests/TestDuplicatePoints.cpp
cs/benchmark/cigma/trunk/src/tests/TestLine.py
cs/benchmark/cigma/trunk/src/tests/TestPathUtils.cpp
cs/benchmark/cigma/trunk/src/tests/TestPointwise.cpp
cs/benchmark/cigma/trunk/src/tests/TestVtkReader.cpp
Modified:
cs/benchmark/cigma/trunk/src/tests/TestLine.cpp
Log:
Test programs
Added: cs/benchmark/cigma/trunk/src/tests/TestDuplicatePoints.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/tests/TestDuplicatePoints.cpp (rev 0)
+++ cs/benchmark/cigma/trunk/src/tests/TestDuplicatePoints.cpp 2008-04-23 16:56:43 UTC (rev 11856)
@@ -0,0 +1,65 @@
+#include <iostream>
+#include <iomanip>
+#include <cmath>
+#include <vector>
+#include <utility>
+#include <cassert>
+
+#include "../PointsReader.h"
+
+using namespace std;
+using namespace cigma;
+
+double norm2(double *x)
+{
+ return x[0]*x[0] + x[1]*x[1] + x[2]*x[2];
+}
+
+double dist2(double *x, double *y)
+{
+ double rho[3] = {x[0]-y[0], x[1]-y[1], x[2]-y[2]};
+ return norm2(rho);
+}
+
+int main()
+{
+ PointsReader pointsReader;
+ //pointsReader.pointsPath = "../tet4_1000m.h5:/projection/points";
+ pointsReader.pointsPath = "./strikeslip_tet4_1000m_t0.vtk:displacements_t0";
+ pointsReader.load_points();
+
+ Points *points = pointsReader.points;
+ if (points == 0)
+ {
+ cerr << "Could not load points!" << endl;
+ exit(1);
+ }
+
+ int npts = points->n_points();
+ int dim = points->n_dim();
+ assert(dim == 3);
+
+ int a,b;
+ const double epsilon2 = 1e-12;
+ for (a = 0; a < npts; a++)
+ {
+ double *x = &(points->data[3*a]);
+ for (b = 0; b < npts; b++)
+ {
+ if (b >= a)
+ {
+ break;
+ }
+
+ double *y = &(points->data[3*b]);
+ double r2 = dist2(x,y);
+ if (r2 < epsilon2)
+ {
+ cout << "(" << a << "," << b << ") ";
+ }
+ }
+ }
+ cout << endl;
+
+ return 0;
+}
Modified: cs/benchmark/cigma/trunk/src/tests/TestLine.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/tests/TestLine.cpp 2008-04-23 16:56:42 UTC (rev 11855)
+++ cs/benchmark/cigma/trunk/src/tests/TestLine.cpp 2008-04-23 16:56:43 UTC (rev 11856)
@@ -22,10 +22,16 @@
int main()
{
- int N = 400;
- double A[3] = { 0, 18000, -12000 };
- double B[3] = { 24000, 18000, -12000 };
+ int N = 1001;
+
+ // line parallel to y-axis
+ //double A[3] = { 0, 12000, -12000 };
+ //double B[3] = { 24000, 12000, -12000 };
+ // line parallel to x-axis
+ double A[3] = { 12010, 0, -12000 };
+ double B[3] = { 12010, 24000, -12000 };
+
int i,j;
double t, dt;
dt = 1.0 / (N - 1);
@@ -64,8 +70,10 @@
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;
@@ -79,7 +87,8 @@
cerr << "Expecting a point field" << endl;
exit(1);
}
- //PointField *pointField = static_cast<PointField*>(fieldReader.field);
+ PointField *pointField = static_cast<PointField*>(fieldReader.field);
+ //pointField->locator = 0; // clear out locator -- brute force search
for (i = 0; i < N; i++)
{
double *station = &x[3*i];
Added: cs/benchmark/cigma/trunk/src/tests/TestLine.py
===================================================================
--- cs/benchmark/cigma/trunk/src/tests/TestLine.py (rev 0)
+++ cs/benchmark/cigma/trunk/src/tests/TestLine.py 2008-04-23 16:56:43 UTC (rev 11856)
@@ -0,0 +1,23 @@
+#!/usr/bin/env python
+
+import numpy
+
+# read array from temporary file
+asdf = numpy.fromfile('asdf', sep=' ')
+asdf = asdf.reshape((1001, 17))
+
+# assign meaning to each column
+n = asdf[:,0]
+t = asdf[:,1]
+x = asdf[:,2:5]
+y0 = asdf[:,5:8]
+y1 = asdf[:,8:11]
+rho = asdf[:,11:14]
+a = asdf[:,14]
+b = asdf[:,15]
+r = asdf[:,16]
+
+# pylab example
+#plot(t,y0,t,y1)
+#legend(('dx','dy','dz','dx2','dy2','dz2'))
+
Added: cs/benchmark/cigma/trunk/src/tests/TestPathUtils.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/tests/TestPathUtils.cpp (rev 0)
+++ cs/benchmark/cigma/trunk/src/tests/TestPathUtils.cpp 2008-04-23 16:56:43 UTC (rev 11856)
@@ -0,0 +1,12 @@
+#include <iostream>
+#include <cassert>
+#include "../PathUtils.h"
+
+int main(int argc, char *argv[])
+{
+
+ filename = argv[1];
+ bool exists = FileExists(filename);
+
+ return 0;
+}
Copied: cs/benchmark/cigma/trunk/src/tests/TestPointwise.cpp (from rev 11855, cs/benchmark/cigma/trunk/src/tests/TestLine.cpp)
===================================================================
--- cs/benchmark/cigma/trunk/src/tests/TestPointwise.cpp (rev 0)
+++ cs/benchmark/cigma/trunk/src/tests/TestPointwise.cpp 2008-04-23 16:56:43 UTC (rev 11856)
@@ -0,0 +1,129 @@
+#include <iostream>
+#include <iomanip>
+#include <cmath>
+#include <cassert>
+
+#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 i,j;
+ int N0 = 1200;
+ int N = 1000;
+
+
+ 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*>(field);
+ Points *points = pointField->points;
+ Points *values = pointField->values;
+
+ int npts = points->n_points();
+ assert(npts > N);
+
+ double *x = points->data;
+
+ 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);
+ }
+
+ for (i = 0; i < N; i++)
+ {
+ for (j = 0; j < 3; j++)
+ {
+ y[1][3*i + j] = values->data[3*(i + N0) + j];
+ }
+ }
+
+ 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];
+ }
+ }
+
+ double epsilon = 0.0;
+ cout << setprecision(12);
+ for (i = 0; i < N; i++)
+ {
+ cout << i << " ";
+ double *X = &x[3*i];
+ double *Y0 = &y[0][3*i];
+ double *Y1 = &y[1][3*i];
+ double *D = &residuals[3*i];
+ double eps = magnitude(D);
+ epsilon += eps * eps;
+
+ print_point(X);
+ print_point(Y0);
+ print_point(Y1);
+ print_point(D);
+
+ cout << magnitude(Y0) << " ";
+ cout << magnitude(Y1) << " ";
+ cout << eps << " ";
+ cout << endl;
+ }
+
+ cerr << setprecision(12);
+ cerr << "epsilon = " << sqrt(epsilon) << endl;
+
+
+ delete [] x;
+ delete [] y[0];
+ delete [] y[1];
+ delete [] residuals;
+
+ return 0;
+}
Added: cs/benchmark/cigma/trunk/src/tests/TestVtkReader.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/tests/TestVtkReader.cpp (rev 0)
+++ cs/benchmark/cigma/trunk/src/tests/TestVtkReader.cpp 2008-04-23 16:56:43 UTC (rev 11856)
@@ -0,0 +1,25 @@
+#include <iostream>
+#include <cassert>
+#include "../VtkReader.h"
+
+using namespace std;
+using namespace cigma;
+
+
+int main(int argc, char *argv[])
+{
+ if (argc == 1)
+ {
+ cerr << "Usage: " << argv[0] << " vtk-file" << endl;
+ return -1;
+ }
+
+ string filename = argv[1];
+ cout << "Reading file '" << filename << "'\n";
+
+ VtkReader reader;
+
+ reader.open(filename.c_str());
+
+ return 0;
+}
More information about the cig-commits
mailing list