[cig-commits] r7498 - in cs/cigma/branches/cigma-0.9/benchmarks: .
disloc3d
luis at geodynamics.org
luis at geodynamics.org
Tue Jun 26 09:41:33 PDT 2007
Author: luis
Date: 2007-06-26 09:41:33 -0700 (Tue, 26 Jun 2007)
New Revision: 7498
Added:
cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/
cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/README
cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/cervelli.c
cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/cervelli.h
cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/dc3d.f
cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/dc3d.h
cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/disloc3d.c
cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/disloc3d.h
cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/h5array.c
cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/h5array.h
cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/okadasoln_ng.c
cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/txtarray.c
cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/txtarray.h
Log:
darcs patch: Added okada solution
Added: cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/README
===================================================================
Added: cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/cervelli.c
===================================================================
--- cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/cervelli.c 2007-06-26 16:40:29 UTC (rev 7497)
+++ cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/cervelli.c 2007-06-26 16:41:33 UTC (rev 7498)
@@ -0,0 +1,247 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "disloc3d.h"
+#include "cervelli.h"
+
+/* function fault_params_to_okada_form
+ *
+ * Based on fault_params_to_okada_form.m MATLAB script
+ * (comments copied verbatim)
+ *
+ * This function takes fault trace, dip and locking depth information
+ * and calculates the anchor coordinates, length, width and strike
+ * of the fault plane as per Okada (1985).
+ *
+ * It may be beneficial in the future to compute the midpoint base
+ * anchor as well for using Okada (1992) Stanford binary.
+ *
+ * We should also allow for buried faults in the future. This can
+ * be achieved by passing a locking depth and a burial depth. The
+ * only output changed would be the width of the fault plane.
+ *
+ * We may need to be more formal about say endpoint1 < endpoint2
+ * ... maybe
+ *
+ * Arguments:
+ * sx1: x coord of fault trace endpoint1
+ * sy1: y coord of fault trace endpoint1
+ * sx2: x coord of fault trace endpoint2
+ * sy2: y coord of fault trace endpoint2
+ * dip: dip of fault plane [radians]
+ * D : fault locking depth
+ * bd : burial depth (top "locking depth")
+ *
+ * Returned variables:
+ * strike: strike of fault plane
+ * L : fault length
+ * W : fault width
+ * ofx : x coord of fault anchor
+ * ofy : y coord of fault anchor
+ * ofxe : x coord of other buried corner
+ * ofye : y coord of other buried corner
+ * tfx : x coord of fault anchor (top relative)
+ * tfy : y coord of fault anchor (top relative)
+ * tfxe : x coord of other buried corner (top relative)
+ * tfye : y coord of other buried corner (top relative)
+ */
+static void
+fault_params_to_okada_form(double sx1, double sy1,
+ double sx2, double sy2,
+ double dip, double D, double bd,
+ double *strike,
+ double *L, double *W,
+ double *ofx, double *ofy,
+ double *ofxe, double *ofye,
+ double *tfx, double *tfy,
+ double *tfxe, double *tfye)
+{
+ //
+ // Do some basic checks
+ //
+ if (D < bd)
+ {
+ fprintf(stderr, "Burial depth is greater than locking depth!\n");
+ fprintf(stderr, "Fault width will be negative!\n");
+ }
+ if (D == bd)
+ {
+ fprintf(stderr, "Burial depth is equal to locking depth!\n");
+ fprintf(stderr, "Fault width will be zero!\n");
+ }
+
+ //
+ // Calculate needed parameters for Okada input
+ //
+ *strike = atan2(sy1 - sy2, sx1 - sx2) + M_PI;
+ *L = sqrt((sx2 - sx1)*(sx2 - sx1) + (sy2 - sy1)*(sy2 - sy1));
+ *W = (D - bd) / sin(dip);
+
+ //
+ // For older versions without a burial depth option
+ //
+ // *W = D/sin(dip)
+ //
+
+ //
+ // Calculate fault segment anchor and other buried point
+ //
+ *ofx = sx1 + D / tan(dip) + sin(*strike);
+ *ofy = sy1 - D / tan(dip) + cos(*strike);
+ *ofxe = sx2 + D / tan(dip) + sin(*strike);
+ *ofye = sy2 - D / tan(dip) + cos(*strike);
+
+ //
+ // Calculate fault segment anchor and other buried point (top relative)
+ //
+ *tfx = sx1 + bd / tan(dip) + sin(*strike);
+ *tfy = sy1 - bd / tan(dip) + cos(*strike);
+ *tfxe = sx2 + bd / tan(dip) + sin(*strike);
+ *tfye = sy2 - bd / tan(dip) + cos(*strike);
+}
+
+
+
+/* function getM
+ * Convert fault parameters into form needed for Cervelli calculations.
+ *
+ * Inputs: fi fault info
+ * nmod number of fault patches
+ * mu Poisson's ratio
+ *
+ * Outputs: M matrix with fault stuff
+ * (see below for what is on the rows)
+ * (one row per fault patch)
+ */
+static void
+getM(double *fi, int nmod, double mu, double *M)
+{
+ int j;
+ double *fault;
+ double *model;
+
+ double fx1,fy1;
+ double fx2,fy2;
+ double dip;
+ double D;
+ double Pr;
+ double bd; // burial depth
+ double ss;
+ double ds;
+ double ts;
+
+ double strike;
+ double L, W;
+ double ofx, ofy;
+ double ofxe, ofye;
+ double tfx, tfy;
+ double tfxe, tfye;
+
+ double length;
+ double width;
+ double depth;
+ double str;
+ double east;
+ double north;
+
+ for (j = 0; j < nmod; j++)
+ {
+ fault = &fi[11*j];
+ model = &M[10*j];
+
+ fx1 = fault[0];
+ fy1 = fault[1];
+ fx2 = fault[2];
+ fy2 = fault[3];
+ dip = fault[4]; // in degs here
+ D = fault[5];
+ Pr = fault[6];
+ bd = fault[7];
+ ss = fault[8];
+ ds = fault[9];
+ ts = fault[10];
+
+ fault_params_to_okada_form(
+ fx1, fy1, fx2, fy2, dip*M_PI/180.0, D, bd,
+ &strike, &L, &W, &ofx, &ofy, &ofxe, &ofye, &tfx, &tfy, &tfxe, &tfye);
+
+ length = L;
+ width = W;
+
+ if (dip < 0)
+ {
+ east = (L/2) * cos(strike) + tfx; // east component of top center
+ north = (L/2) * sin(strike) + tfy; // north " " " "
+ depth = bd;
+ width = abs(width);
+ ss = -1 * ss;
+ ds = -1 * ds;
+ ts = -1 * ts;
+ }
+ else
+ {
+ east = (L/2) * cos(strike) + ofx; // east component of bottom center
+ north = (L/2) * sin(strike) + ofy;
+ depth = D;
+ }
+
+ // Cervelli takes strike in degs (strike comes out of
+ // fault_params_to_okada_form in rads)
+ str = (180.0/M_PI) * ((M_PI/2) - strike);
+
+ model[0] = length;
+ model[1] = width;
+ model[2] = depth;
+ model[3] = dip;
+ model[4] = str;
+ model[5] = east;
+ model[6] = north;
+ model[7] = ss;
+ model[8] = ds;
+ model[9] = ts;
+ }
+
+}
+
+
+
+/* function calc_disp_cervelli
+ *
+ * Calculate displacement's via Cervelli's Okada (1992) code
+ *
+ * Inputs: fault_info matrix with fault stuff
+ * sx station x coord
+ * sy station y coord
+ * sz station z coord
+ *
+ * Outputs: dispx disp in x direction
+ * dispy disp in y direction
+ * dispz disp in z direction
+ * D matrix with derivatives
+ * S Stresses (not used)
+ *
+ * TODO: update this comment
+ */
+void calc_disp_cervelli(double mu, double nu,
+ double *models, double *fault_info, int nfaults,
+ double *stations, int nstations,
+ double *disp, double *deriv, double *stress,
+ double *flag, double *flag2)
+{
+ //
+ // get the "model matrix" with
+ // fault patch info in the form for cervelli calcs
+ //
+ getM(fault_info, nfaults, mu, models);
+
+ //
+ // call the disloc3d wrapper
+ //
+ disloc3d(models, nfaults,
+ stations, nstations,
+ mu, nu,
+ disp, deriv, stress,
+ flag, flag2);
+
+}
+
Added: cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/cervelli.h
===================================================================
--- cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/cervelli.h 2007-06-26 16:40:29 UTC (rev 7497)
+++ cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/cervelli.h 2007-06-26 16:41:33 UTC (rev 7498)
@@ -0,0 +1,11 @@
+#ifndef __CERVELLI_H__
+#define __CERVELLI_H__
+
+
+void calc_disp_cervelli(double mu, double nu,
+ double *models, double *fault_info, int nfaults,
+ double *stations, int nstations,
+ double *disp, double *deriv, double *stress,
+ double *flag, double *flag2);
+
+#endif
Added: cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/dc3d.f
===================================================================
--- cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/dc3d.f 2007-06-26 16:40:29 UTC (rev 7497)
+++ cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/dc3d.f 2007-06-26 16:41:33 UTC (rev 7498)
@@ -0,0 +1,612 @@
+ SUBROUTINE DC3D(ALPHA,X,Y,Z,DEPTH,DIP, 04610000
+ * AL1,AL2,AW1,AW2,DISL1,DISL2,DISL3, 04620002
+ * UX,UY,UZ,UXX,UYX,UZX,UXY,UYY,UZY,UXZ,UYZ,UZZ,IRET) 04630002
+ IMPLICIT REAL*8 (A-H,O-Z) 04640000
+ REAL*8 ALPHA,X,Y,Z,DEPTH,DIP,AL1,AL2,AW1,AW2,DISL1,DISL2,DISL3, 04650000
+ * UX,UY,UZ,UXX,UYX,UZX,UXY,UYY,UZY,UXZ,UYZ,UZZ,IRET 04660000
+C 04670000
+C******************************************************************** 04680000
+C***** ***** 04690000
+C***** DISPLACEMENT AND STRAIN AT DEPTH ***** 04700000
+C***** DUE TO BURIED FINITE FAULT IN A SEMIINFINITE MEDIUM ***** 04710000
+C***** CODED BY Y.OKADA ... SEP 1991 ***** 04720002
+C***** REVISED Y.OKADA ... NOV 1991 ***** 04730002
+C***** ***** 04740000
+C******************************************************************** 04750000
+C 04760000
+C***** INPUT 04770000
+C***** ALPHA : MEDIUM CONSTANT (LAMBDA+MYU)/(LAMBDA+2*MYU) 04780000
+C***** X,Y,Z : COORDINATE OF OBSERVING POINT 04790000
+C***** DEPTH : SOURCE DEPTH 04800000
+C***** DIP : DIP-ANGLE (DEGREE) 04810000
+C***** AL1,AL2 : FAULT LENGTH (-STRIKE,+STRIKE) 04820000
+C***** AW1,AW2 : FAULT WIDTH ( DOWNDIP, UPDIP) 04830000
+C***** DISL1-DISL3 : STRIKE-, DIP-, TENSILE-DISLOCATIONS 04840000
+C 04850000
+C***** OUTPUT 04860000
+C***** UX, UY, UZ : DISPLACEMENT ( UNIT=(UNIT OF DISL) 04870000
+C***** UXX,UYX,UZX : X-DERIVATIVE ( UNIT=(UNIT OF DISL) / 04880000
+C***** UXY,UYY,UZY : Y-DERIVATIVE (UNIT OF X,Y,Z,DEPTH,AL,AW) )04890000
+C***** UXZ,UYZ,UZZ : Z-DERIVATIVE 04900000
+C***** IRET : RETURN CODE ( =0....NORMAL, =1....SINGULAR ) 04910002
+C 04920000
+ COMMON /C0/DUMMY(5),SD,CD ,Dummy2(5)
+ COMMON /C2/XI2,ET2,Q2,R ,Dummy3(20)
+ DIMENSION U(12),DU(12),DUA(12),DUB(12),DUC(12) 04950000
+ DATA F0/0.D0/ 04960000
+C----- 04970000
+ IF(Z.GT.0.) WRITE(6,'('' ** POSITIVE Z WAS GIVEN IN SUB-DC3D'')') 04980000
+ DO 111 I=1,12 04990000
+ U (I)=F0 05000000
+ DUA(I)=F0 05010000
+ DUB(I)=F0 05020000
+ DUC(I)=F0 05030000
+ 111 CONTINUE 05040000
+ AALPHA=ALPHA 05050000
+ DDIP=DIP 05060000
+ CALL DCCON0(AALPHA,DDIP) 05070000
+C====================================== 05080000
+C===== REAL-SOURCE CONTRIBUTION ===== 05090000
+C====================================== 05100000
+ D=DEPTH+Z 05110000
+ P=Y*CD+D*SD 05120000
+ Q=Y*SD-D*CD 05130000
+ JXI=0 05140000
+ JET=0 05150000
+ IF((X+AL1)*(X-AL2).LE.0.) JXI=1 05160000
+ IF((P+AW1)*(P-AW2).LE.0.) JET=1 05170000
+ DD1=DISL1 05180000
+ DD2=DISL2 05190000
+ DD3=DISL3 05200000
+C----- 05210000
+ DO 223 K=1,2 05220000
+ IF(K.EQ.1) ET=P+AW1 05230000
+ IF(K.EQ.2) ET=P-AW2 05240000
+ DO 222 J=1,2 05250000
+ IF(J.EQ.1) XI=X+AL1 05260000
+ IF(J.EQ.2) XI=X-AL2 05270000
+ CALL DCCON2(XI,ET,Q,SD,CD) 05280002
+ IF(JXI.EQ.1 .AND. Q.EQ.F0 .AND. ET.EQ.F0) GO TO 99 05290002
+ IF(JET.EQ.1 .AND. Q.EQ.F0 .AND. XI.EQ.F0) GO TO 99 05300002
+ CALL UA(XI,ET,Q,DD1,DD2,DD3,DUA) 05310000
+C----- 05320000
+ DO 220 I=1,10,3 05330000
+ DU(I) =-DUA(I) 05340000
+ DU(I+1)=-DUA(I+1)*CD+DUA(I+2)*SD 05350000
+ DU(I+2)=-DUA(I+1)*SD-DUA(I+2)*CD 05360000
+ IF(I.LT.10) GO TO 220 05370000
+ DU(I) =-DU(I) 05380000
+ DU(I+1)=-DU(I+1) 05390000
+ DU(I+2)=-DU(I+2) 05400000
+ 220 CONTINUE 05410000
+ DO 221 I=1,12 05420000
+ IF(J+K.NE.3) U(I)=U(I)+DU(I) 05430000
+ IF(J+K.EQ.3) U(I)=U(I)-DU(I) 05440000
+ 221 CONTINUE 05450000
+C----- 05460000
+ 222 CONTINUE 05470000
+ 223 CONTINUE 05480000
+C======================================= 05490000
+C===== IMAGE-SOURCE CONTRIBUTION ===== 05500000
+C======================================= 05510000
+ ZZ=Z 05520000
+ D=DEPTH-Z 05530000
+ P=Y*CD+D*SD 05540000
+ Q=Y*SD-D*CD 05550000
+ JET=0 05560000
+ IF((P+AW1)*(P-AW2).LE.0.) JET=1 05570000
+C----- 05580000
+ DO 334 K=1,2 05590000
+ IF(K.EQ.1) ET=P+AW1 05600000
+ IF(K.EQ.2) ET=P-AW2 05610000
+ DO 333 J=1,2 05620000
+ IF(J.EQ.1) XI=X+AL1 05630000
+ IF(J.EQ.2) XI=X-AL2 05640000
+ CALL DCCON2(XI,ET,Q,SD,CD) 05650002
+ CALL UA(XI,ET,Q,DD1,DD2,DD3,DUA) 05660000
+ CALL UB(XI,ET,Q,DD1,DD2,DD3,DUB) 05670000
+ CALL UC(XI,ET,Q,ZZ,DD1,DD2,DD3,DUC) 05680000
+C----- 05690000
+ DO 330 I=1,10,3 05700000
+ DU(I)=DUA(I)+DUB(I)+Z*DUC(I) 05710000
+ DU(I+1)=(DUA(I+1)+DUB(I+1)+Z*DUC(I+1))*CD 05720000
+ * -(DUA(I+2)+DUB(I+2)+Z*DUC(I+2))*SD 05730000
+ DU(I+2)=(DUA(I+1)+DUB(I+1)-Z*DUC(I+1))*SD 05740000
+ * +(DUA(I+2)+DUB(I+2)-Z*DUC(I+2))*CD 05750000
+ IF(I.LT.10) GO TO 330 05760000
+ DU(10)=DU(10)+DUC(1) 05770000
+ DU(11)=DU(11)+DUC(2)*CD-DUC(3)*SD 05780000
+ DU(12)=DU(12)-DUC(2)*SD-DUC(3)*CD 05790000
+ 330 CONTINUE 05800000
+ DO 331 I=1,12 05810000
+ IF(J+K.NE.3) U(I)=U(I)+DU(I) 05820000
+ IF(J+K.EQ.3) U(I)=U(I)-DU(I) 05830000
+ 331 CONTINUE 05840000
+C----- 05850000
+ 333 CONTINUE 05860000
+ 334 CONTINUE 05870000
+C===== 05880000
+ UX=U(1) 05890000
+ UY=U(2) 05900000
+ UZ=U(3) 05910000
+ UXX=U(4) 05920000
+ UYX=U(5) 05930000
+ UZX=U(6) 05940000
+ UXY=U(7) 05950000
+ UYY=U(8) 05960000
+ UZY=U(9) 05970000
+ UXZ=U(10) 05980000
+ UYZ=U(11) 05990000
+ UZZ=U(12) 06000000
+ IRET=0 06010002
+ RETURN 06020000
+C======================================= 06030000
+C===== IN CASE OF SINGULAR (R=0) ===== 06040000
+C======================================= 06050000
+ 99 UX=F0 06060000
+ UY=F0 06070000
+ UZ=F0 06080000
+ UXX=F0 06090000
+ UYX=F0 06100000
+ UZX=F0 06110000
+ UXY=F0 06120000
+ UYY=F0 06130000
+ UZY=F0 06140000
+ UXZ=F0 06150000
+ UYZ=F0 06160000
+ UZZ=F0 06170000
+ IRET=1 06180002
+ RETURN 06190000
+ END 06200000
+
+
+
+ SUBROUTINE DCCON0(ALPHA,DIP) 09300000
+ IMPLICIT REAL*8 (A-H,O-Z) 09310000
+C 09320000
+C******************************************************************* 09330000
+C***** CALCULATE MEDIUM CONSTANTS AND FAULT-DIP CONSTANTS ***** 09340000
+C******************************************************************* 09350000
+C 09360000
+C***** INPUT 09370000
+C***** ALPHA : MEDIUM CONSTANT (LAMBDA+MYU)/(LAMBDA+2*MYU) 09380000
+C***** DIP : DIP-ANGLE (DEGREE) 09390000
+C### CAUTION ### IF COS(DIP) IS SUFFICIENTLY SMALL, IT IS SET TO ZERO 09400000
+C 09410000
+ COMMON /C0/ALP1,ALP2,ALP3,ALP4,ALP5,SD,CD,SDSD,CDCD,SDCD,S2D,C2D 09420000
+c everything in CO is used in this routine
+
+ DATA F0,F1,F2,PI2/0.D0,1.D0,2.D0,6.283185307179586D0/ 09430000
+ DATA EPS/1.D-6/ 09440000
+C----- 09450000
+ ALP1=(F1-ALPHA)/F2 09460000
+ ALP2= ALPHA/F2 09470000
+ ALP3=(F1-ALPHA)/ALPHA 09480000
+ ALP4= F1-ALPHA 09490000
+ ALP5= ALPHA 09500000
+C----- 09510000
+ P18=PI2/360.D0 09520000
+ SD=DSIN(DIP*P18) 09530000
+ CD=DCOS(DIP*P18) 09540000
+ IF(DABS(CD).LT.EPS) THEN 09550000
+ CD=F0 09560000
+ IF(SD.GT.F0) SD= F1 09570000
+ IF(SD.LT.F0) SD=-F1 09580000
+ ENDIF 09590000
+ SDSD=SD*SD 09600000
+ CDCD=CD*CD 09610000
+ SDCD=SD*CD 09620000
+ S2D=F2*SDCD 09630000
+ C2D=CDCD-SDSD 09640000
+ RETURN 09650000
+ END 09660000
+
+
+ SUBROUTINE UA(XI,ET,Q,DISL1,DISL2,DISL3,U) 06210000
+ IMPLICIT REAL*8 (A-H,O-Z) 06220000
+ DIMENSION U(12),DU(12) 06230000
+C 06240000
+C******************************************************************** 06250000
+C***** DISPLACEMENT AND STRAIN AT DEPTH (PART-A) ***** 06260000
+C***** DUE TO BURIED FINITE FAULT IN A SEMIINFINITE MEDIUM ***** 06270000
+C******************************************************************** 06280000
+C 06290000
+C***** INPUT 06300000
+C***** XI,ET,Q : STATION COORDINATES IN FAULT SYSTEM 06310000
+C***** DISL1-DISL3 : STRIKE-, DIP-, TENSILE-DISLOCATIONS 06320000
+C***** OUTPUT 06330000
+C***** U(12) : DISPLACEMENT AND THEIR DERIVATIVES 06340000
+C 06350000
+ COMMON /C0/ALP1,ALP2,ALP3,ALP4,ALP5,SD,CD,SDSD,CDCD,SDCD,S2D,C2D 06360000
+ COMMON /C2/XI2,ET2,Q2,R,R2,R3,R5,Y,D,TT,ALX,ALE,X11,Y11,X32,Y32, 06370000
+ * EY,EZ,FY,FZ,GY,GZ,HY,HZ 06380000
+ DATA F0,F2,PI2/0.D0,2.D0,6.283185307179586D0/ 06390000
+C----- 06400000
+ DO 111 I=1,12 06410000
+ 111 U(I)=F0 06420000
+ XY=XI*Y11 06430000
+ QX=Q *X11 06440000
+ QY=Q *Y11 06450000
+C====================================== 06460000
+C===== STRIKE-SLIP CONTRIBUTION ===== 06470000
+C====================================== 06480000
+ IF(DISL1.NE.F0) THEN 06490000
+ DU( 1)= TT/F2 +ALP2*XI*QY 06500000
+ DU( 2)= ALP2*Q/R 06510000
+ DU( 3)= ALP1*ALE -ALP2*Q*QY 06520000
+ DU( 4)=-ALP1*QY -ALP2*XI2*Q*Y32 06530000
+ DU( 5)= -ALP2*XI*Q/R3 06540000
+ DU( 6)= ALP1*XY +ALP2*XI*Q2*Y32 06550000
+ DU( 7)= ALP1*XY*SD +ALP2*XI*FY+D/F2*X11 06560000
+ DU( 8)= ALP2*EY 06570000
+ DU( 9)= ALP1*(CD/R+QY*SD) -ALP2*Q*FY 06580000
+ DU(10)= ALP1*XY*CD +ALP2*XI*FZ+Y/F2*X11 06590000
+ DU(11)= ALP2*EZ 06600000
+ DU(12)=-ALP1*(SD/R-QY*CD) -ALP2*Q*FZ 06610000
+ DO 222 I=1,12 06620000
+ 222 U(I)=U(I)+DISL1/PI2*DU(I) 06630000
+ ENDIF 06640000
+C====================================== 06650000
+C===== DIP-SLIP CONTRIBUTION ===== 06660000
+C====================================== 06670000
+ IF(DISL2.NE.F0) THEN 06680000
+ DU( 1)= ALP2*Q/R 06690000
+ DU( 2)= TT/F2 +ALP2*ET*QX 06700000
+ DU( 3)= ALP1*ALX -ALP2*Q*QX 06710000
+ DU( 4)= -ALP2*XI*Q/R3 06720000
+ DU( 5)= -QY/F2 -ALP2*ET*Q/R3 06730000
+ DU( 6)= ALP1/R +ALP2*Q2/R3 06740000
+ DU( 7)= ALP2*EY 06750000
+ DU( 8)= ALP1*D*X11+XY/F2*SD +ALP2*ET*GY 06760000
+ DU( 9)= ALP1*Y*X11 -ALP2*Q*GY 06770000
+ DU(10)= ALP2*EZ 06780000
+ DU(11)= ALP1*Y*X11+XY/F2*CD +ALP2*ET*GZ 06790000
+ DU(12)=-ALP1*D*X11 -ALP2*Q*GZ 06800000
+ DO 333 I=1,12 06810000
+ 333 U(I)=U(I)+DISL2/PI2*DU(I) 06820000
+ ENDIF 06830000
+C======================================== 06840000
+C===== TENSILE-FAULT CONTRIBUTION ===== 06850000
+C======================================== 06860000
+ IF(DISL3.NE.F0) THEN 06870000
+ DU( 1)=-ALP1*ALE -ALP2*Q*QY 06880000
+ DU( 2)=-ALP1*ALX -ALP2*Q*QX 06890000
+ DU( 3)= TT/F2 -ALP2*(ET*QX+XI*QY) 06900000
+ DU( 4)=-ALP1*XY +ALP2*XI*Q2*Y32 06910000
+ DU( 5)=-ALP1/R +ALP2*Q2/R3 06920000
+ DU( 6)=-ALP1*QY -ALP2*Q*Q2*Y32 06930000
+ DU( 7)=-ALP1*(CD/R+QY*SD) -ALP2*Q*FY 06940000
+ DU( 8)=-ALP1*Y*X11 -ALP2*Q*GY 06950000
+ DU( 9)= ALP1*(D*X11+XY*SD) +ALP2*Q*HY 06960000
+ DU(10)= ALP1*(SD/R-QY*CD) -ALP2*Q*FZ 06970000
+ DU(11)= ALP1*D*X11 -ALP2*Q*GZ 06980000
+ DU(12)= ALP1*(Y*X11+XY*CD) +ALP2*Q*HZ 06990000
+ DO 444 I=1,12 07000000
+ 444 U(I)=U(I)+DISL3/PI2*DU(I) 07010000
+ ENDIF 07020000
+ RETURN 07030000
+ END 07040000
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ SUBROUTINE UB(XI,ET,Q,DISL1,DISL2,DISL3,U) 07050000
+ IMPLICIT REAL*8 (A-H,O-Z) 07060000
+ DIMENSION U(12),DU(12) 07070000
+C 07080000
+C******************************************************************** 07090000
+C***** DISPLACEMENT AND STRAIN AT DEPTH (PART-B) ***** 07100000
+C***** DUE TO BURIED FINITE FAULT IN A SEMIINFINITE MEDIUM ***** 07110000
+C******************************************************************** 07120000
+C 07130000
+C***** INPUT 07140000
+C***** XI,ET,Q : STATION COORDINATES IN FAULT SYSTEM 07150000
+C***** DISL1-DISL3 : STRIKE-, DIP-, TENSILE-DISLOCATIONS 07160000
+C***** OUTPUT 07170000
+C***** U(12) : DISPLACEMENT AND THEIR DERIVATIVES 07180000
+C 07190000
+ COMMON /C0/ALP1,ALP2,ALP3,ALP4,ALP5,SD,CD,SDSD,CDCD,SDCD,S2D,C2D 07200000
+ COMMON /C2/XI2,ET2,Q2,R,R2,R3,R5,Y,D,TT,ALX,ALE,X11,Y11,X32,Y32, 07210000
+ * EY,EZ,FY,FZ,GY,GZ,HY,HZ 07220000
+ DATA F0,F1,F2,PI2/0.D0,1.D0,2.D0,6.283185307179586D0/ 07230000
+C----- 07240000
+ RD=R+D 07250000
+ D11=F1/(R*RD) 07260000
+ AJ2=XI*Y/RD*D11 07270000
+ AJ5=-(D+Y*Y/RD)*D11 07280000
+ IF(CD.NE.F0) THEN 07290000
+ IF(XI.EQ.F0) THEN 07300000
+ AI4=F0 07310000
+ ELSE 07320000
+ X=DSQRT(XI2+Q2) 07330000
+ AI4=F1/CDCD*( XI/RD*SDCD 07340000
+ * +F2*DATAN((ET*(X+Q*CD)+X*(R+X)*SD)/(XI*(R+X)*CD)) ) 07350000
+ ENDIF 07360000
+ AI3=(Y*CD/RD-ALE+SD*DLOG(RD))/CDCD 07370000
+ AK1=XI*(D11-Y11*SD)/CD 07380000
+ AK3=(Q*Y11-Y*D11)/CD 07390000
+ AJ3=(AK1-AJ2*SD)/CD 07400000
+ AJ6=(AK3-AJ5*SD)/CD 07410000
+ ELSE 07420000
+ RD2=RD*RD 07430000
+ AI3=(ET/RD+Y*Q/RD2-ALE)/F2 07440000
+ AI4=XI*Y/RD2/F2 07450000
+ AK1=XI*Q/RD*D11 07460000
+ AK3=SD/RD*(XI2*D11-F1) 07470000
+ AJ3=-XI/RD2*(Q2*D11-F1/F2) 07480000
+ AJ6=-Y/RD2*(XI2*D11-F1/F2) 07490000
+ ENDIF 07500000
+C----- 07510000
+ XY=XI*Y11 07520000
+ AI1=-XI/RD*CD-AI4*SD 07530000
+ AI2= DLOG(RD)+AI3*SD 07540000
+ AK2= F1/R+AK3*SD 07550000
+ AK4= XY*CD-AK1*SD 07560000
+ AJ1= AJ5*CD-AJ6*SD 07570000
+ AJ4=-XY-AJ2*CD+AJ3*SD 07580000
+C===== 07590000
+ DO 111 I=1,12 07600000
+ 111 U(I)=F0 07610000
+ QX=Q*X11 07620000
+ QY=Q*Y11 07630000
+C====================================== 07640000
+C===== STRIKE-SLIP CONTRIBUTION ===== 07650000
+C====================================== 07660000
+ IF(DISL1.NE.F0) THEN 07670000
+ DU( 1)=-XI*QY-TT -ALP3*AI1*SD 07680000
+ DU( 2)=-Q/R +ALP3*Y/RD*SD 07690000
+ DU( 3)= Q*QY -ALP3*AI2*SD 07700000
+ DU( 4)= XI2*Q*Y32 -ALP3*AJ1*SD 07710000
+ DU( 5)= XI*Q/R3 -ALP3*AJ2*SD 07720000
+ DU( 6)=-XI*Q2*Y32 -ALP3*AJ3*SD 07730000
+ DU( 7)=-XI*FY-D*X11 +ALP3*(XY+AJ4)*SD 07740000
+ DU( 8)=-EY +ALP3*(F1/R+AJ5)*SD 07750000
+ DU( 9)= Q*FY -ALP3*(QY-AJ6)*SD 07760000
+ DU(10)=-XI*FZ-Y*X11 +ALP3*AK1*SD 07770000
+ DU(11)=-EZ +ALP3*Y*D11*SD 07780000
+ DU(12)= Q*FZ +ALP3*AK2*SD 07790000
+ DO 222 I=1,12 07800000
+ 222 U(I)=U(I)+DISL1/PI2*DU(I) 07810000
+ ENDIF 07820000
+C====================================== 07830000
+C===== DIP-SLIP CONTRIBUTION ===== 07840000
+C====================================== 07850000
+ IF(DISL2.NE.F0) THEN 07860000
+ DU( 1)=-Q/R +ALP3*AI3*SDCD 07870000
+ DU( 2)=-ET*QX-TT -ALP3*XI/RD*SDCD 07880000
+ DU( 3)= Q*QX +ALP3*AI4*SDCD 07890000
+ DU( 4)= XI*Q/R3 +ALP3*AJ4*SDCD 07900000
+ DU( 5)= ET*Q/R3+QY +ALP3*AJ5*SDCD 07910000
+ DU( 6)=-Q2/R3 +ALP3*AJ6*SDCD 07920000
+ DU( 7)=-EY +ALP3*AJ1*SDCD 07930000
+ DU( 8)=-ET*GY-XY*SD +ALP3*AJ2*SDCD 07940000
+ DU( 9)= Q*GY +ALP3*AJ3*SDCD 07950000
+ DU(10)=-EZ -ALP3*AK3*SDCD 07960000
+ DU(11)=-ET*GZ-XY*CD -ALP3*XI*D11*SDCD 07970000
+ DU(12)= Q*GZ -ALP3*AK4*SDCD 07980000
+ DO 333 I=1,12 07990000
+ 333 U(I)=U(I)+DISL2/PI2*DU(I) 08000000
+ ENDIF 08010000
+C======================================== 08020000
+C===== TENSILE-FAULT CONTRIBUTION ===== 08030000
+C======================================== 08040000
+ IF(DISL3.NE.F0) THEN 08050000
+ DU( 1)= Q*QY -ALP3*AI3*SDSD 08060000
+ DU( 2)= Q*QX +ALP3*XI/RD*SDSD 08070000
+ DU( 3)= ET*QX+XI*QY-TT -ALP3*AI4*SDSD 08080000
+ DU( 4)=-XI*Q2*Y32 -ALP3*AJ4*SDSD 08090000
+ DU( 5)=-Q2/R3 -ALP3*AJ5*SDSD 08100000
+ DU( 6)= Q*Q2*Y32 -ALP3*AJ6*SDSD 08110000
+ DU( 7)= Q*FY -ALP3*AJ1*SDSD 08120000
+ DU( 8)= Q*GY -ALP3*AJ2*SDSD 08130000
+ DU( 9)=-Q*HY -ALP3*AJ3*SDSD 08140000
+ DU(10)= Q*FZ +ALP3*AK3*SDSD 08150000
+ DU(11)= Q*GZ +ALP3*XI*D11*SDSD 08160000
+ DU(12)=-Q*HZ +ALP3*AK4*SDSD 08170000
+ DO 444 I=1,12 08180000
+ 444 U(I)=U(I)+DISL3/PI2*DU(I) 08190000
+ ENDIF 08200000
+ RETURN 08210000
+ END 08220000
+
+
+
+
+
+ SUBROUTINE UC(XI,ET,Q,Z,DISL1,DISL2,DISL3,U) 08230000
+ IMPLICIT REAL*8 (A-H,O-Z) 08240000
+ DIMENSION U(12),DU(12) 08250000
+C 08260000
+C******************************************************************** 08270000
+C***** DISPLACEMENT AND STRAIN AT DEPTH (PART-C) ***** 08280000
+C***** DUE TO BURIED FINITE FAULT IN A SEMIINFINITE MEDIUM ***** 08290000
+C******************************************************************** 08300000
+C 08310000
+C***** INPUT 08320000
+C***** XI,ET,Q,Z : STATION COORDINATES IN FAULT SYSTEM 08330000
+C***** DISL1-DISL3 : STRIKE-, DIP-, TENSILE-DISLOCATIONS 08340000
+C***** OUTPUT 08350000
+C***** U(12) : DISPLACEMENT AND THEIR DERIVATIVES 08360000
+C 08370000
+ COMMON /C0/ALP1,ALP2,ALP3,ALP4,ALP5,SD,CD,SDSD,CDCD,SDCD,S2D,C2D 08380000
+ COMMON /C2/XI2,ET2,Q2,R,R2,R3,R5,Y,D,TT,ALX,ALE,X11,Y11,X32,Y32, 08390000
+ * EY,EZ,FY,FZ,GY,GZ,HY,HZ 08400000
+ DATA F0,F1,F2,F3,PI2/0.D0,1.D0,2.D0,3.D0,6.283185307179586D0/ 08410000
+C----- 08420000
+ C=D+Z 08430000
+ X53=(8.D0*R2+9.D0*R*XI+F3*XI2)*X11*X11*X11/R2 08440000
+ Y53=(8.D0*R2+9.D0*R*ET+F3*ET2)*Y11*Y11*Y11/R2 08450000
+ H=Q*CD-Z 08460000
+ Z32=SD/R3-H*Y32 08470000
+ Z53=F3*SD/R5-H*Y53 08480000
+ Y0=Y11-XI2*Y32 08490000
+ Z0=Z32-XI2*Z53 08500000
+ PPY=CD/R3+Q*Y32*SD 08510000
+ PPZ=SD/R3-Q*Y32*CD 08520000
+ QQ=Z*Y32+Z32+Z0 08530000
+ QQY=F3*C*D/R5-QQ*SD 08540000
+ QQZ=F3*C*Y/R5-QQ*CD+Q*Y32 08550000
+ XY=XI*Y11 08560000
+ QX=Q*X11 08570000
+ QY=Q*Y11 08580000
+ QR=F3*Q/R5 08590000
+ CQX=C*Q*X53 08600000
+ CDR=(C+D)/R3 08610000
+ YY0=Y/R3-Y0*CD 08620000
+C===== 08630000
+ DO 111 I=1,12 08640000
+ 111 U(I)=F0 08650000
+C====================================== 08660000
+C===== STRIKE-SLIP CONTRIBUTION ===== 08670000
+C====================================== 08680000
+ IF(DISL1.NE.F0) THEN 08690000
+ DU( 1)= ALP4*XY*CD -ALP5*XI*Q*Z32 08700000
+ DU( 2)= ALP4*(CD/R+F2*QY*SD) -ALP5*C*Q/R3 08710000
+ DU( 3)= ALP4*QY*CD -ALP5*(C*ET/R3-Z*Y11+XI2*Z32) 08720000
+ DU( 4)= ALP4*Y0*CD -ALP5*Q*Z0 08730000
+ DU( 5)=-ALP4*XI*(CD/R3+F2*Q*Y32*SD) +ALP5*C*XI*QR 08740000
+ DU( 6)=-ALP4*XI*Q*Y32*CD +ALP5*XI*(F3*C*ET/R5-QQ) 08750000
+ DU( 7)=-ALP4*XI*PPY*CD -ALP5*XI*QQY 08760000
+ DU( 8)= ALP4*F2*(D/R3-Y0*SD)*SD-Y/R3*CD 08770000
+ * -ALP5*(CDR*SD-ET/R3-C*Y*QR) 08780000
+ DU( 9)=-ALP4*Q/R3+YY0*SD +ALP5*(CDR*CD+C*D*QR-(Y0*CD+Q*Z0)*SD) 08790000
+ DU(10)= ALP4*XI*PPZ*CD -ALP5*XI*QQZ 08800000
+ DU(11)= ALP4*F2*(Y/R3-Y0*CD)*SD+D/R3*CD -ALP5*(CDR*CD+C*D*QR) 08810000
+ DU(12)= YY0*CD -ALP5*(CDR*SD-C*Y*QR-Y0*SDSD+Q*Z0*CD) 08820000
+ DO 222 I=1,12 08830000
+ 222 U(I)=U(I)+DISL1/PI2*DU(I) 08840000
+ ENDIF 08850000
+C====================================== 08860000
+C===== DIP-SLIP CONTRIBUTION ===== 08870000
+C====================================== 08880000
+ IF(DISL2.NE.F0) THEN 08890000
+ DU( 1)= ALP4*CD/R -QY*SD -ALP5*C*Q/R3 08900000
+ DU( 2)= ALP4*Y*X11 -ALP5*C*ET*Q*X32 08910000
+ DU( 3)= -D*X11-XY*SD -ALP5*C*(X11-Q2*X32) 08920000
+ DU( 4)=-ALP4*XI/R3*CD +ALP5*C*XI*QR +XI*Q*Y32*SD 08930000
+ DU( 5)=-ALP4*Y/R3 +ALP5*C*ET*QR 08940000
+ DU( 6)= D/R3-Y0*SD +ALP5*C/R3*(F1-F3*Q2/R2) 08950000
+ DU( 7)=-ALP4*ET/R3+Y0*SDSD -ALP5*(CDR*SD-C*Y*QR) 08960000
+ DU( 8)= ALP4*(X11-Y*Y*X32) -ALP5*C*((D+F2*Q*CD)*X32-Y*ET*Q*X53) 08970000
+ DU( 9)= XI*PPY*SD+Y*D*X32 +ALP5*C*((Y+F2*Q*SD)*X32-Y*Q2*X53) 08980000
+ DU(10)= -Q/R3+Y0*SDCD -ALP5*(CDR*CD+C*D*QR) 08990000
+ DU(11)= ALP4*Y*D*X32 -ALP5*C*((Y-F2*Q*SD)*X32+D*ET*Q*X53) 09000000
+ DU(12)=-XI*PPZ*SD+X11-D*D*X32-ALP5*C*((D-F2*Q*CD)*X32-D*Q2*X53) 09010000
+ DO 333 I=1,12 09020000
+ 333 U(I)=U(I)+DISL2/PI2*DU(I) 09030000
+ ENDIF 09040000
+C======================================== 09050000
+C===== TENSILE-FAULT CONTRIBUTION ===== 09060000
+C======================================== 09070000
+ IF(DISL3.NE.F0) THEN 09080000
+ DU( 1)=-ALP4*(SD/R+QY*CD) -ALP5*(Z*Y11-Q2*Z32) 09090000
+ DU( 2)= ALP4*F2*XY*SD+D*X11 -ALP5*C*(X11-Q2*X32) 09100000
+ DU( 3)= ALP4*(Y*X11+XY*CD) +ALP5*Q*(C*ET*X32+XI*Z32) 09110000
+ DU( 4)= ALP4*XI/R3*SD+XI*Q*Y32*CD+ALP5*XI*(F3*C*ET/R5-F2*Z32-Z0)09120000
+ DU( 5)= ALP4*F2*Y0*SD-D/R3 +ALP5*C/R3*(F1-F3*Q2/R2) 09130000
+ DU( 6)=-ALP4*YY0 -ALP5*(C*ET*QR-Q*Z0) 09140000
+ DU( 7)= ALP4*(Q/R3+Y0*SDCD) +ALP5*(Z/R3*CD+C*D*QR-Q*Z0*SD) 09150000
+ DU( 8)=-ALP4*F2*XI*PPY*SD-Y*D*X32 09160000
+ * +ALP5*C*((Y+F2*Q*SD)*X32-Y*Q2*X53) 09170000
+ DU( 9)=-ALP4*(XI*PPY*CD-X11+Y*Y*X32) 09180000
+ * +ALP5*(C*((D+F2*Q*CD)*X32-Y*ET*Q*X53)+XI*QQY) 09190000
+ DU(10)= -ET/R3+Y0*CDCD -ALP5*(Z/R3*SD-C*Y*QR-Y0*SDSD+Q*Z0*CD) 09200000
+ DU(11)= ALP4*F2*XI*PPZ*SD-X11+D*D*X32 09210000
+ * -ALP5*C*((D-F2*Q*CD)*X32-D*Q2*X53) 09220000
+ DU(12)= ALP4*(XI*PPZ*CD+Y*D*X32) 09230000
+ * +ALP5*(C*((Y-F2*Q*SD)*X32+D*ET*Q*X53)+XI*QQZ) 09240000
+ DO 444 I=1,12 09250000
+ 444 U(I)=U(I)+DISL3/PI2*DU(I) 09260000
+ ENDIF 09270000
+ RETURN 09280000
+ END 09290000
+
+
+
+
+
+
+
+ SUBROUTINE DCCON2(XI,ET,Q,SD,CD) 10170002
+ IMPLICIT REAL*8 (A-H,O-Z) 10180000
+C 10190000
+C********************************************************************** 10200000
+C***** CALCULATE STATION GEOMETRY CONSTANTS FOR FINITE SOURCE ***** 10210000
+C********************************************************************** 10220000
+C 10230000
+C***** INPUT 10240000
+C***** XI,ET,Q : STATION COORDINATES IN FAULT SYSTEM 10250000
+C***** SD,CD : SIN, COS OF DIP-ANGLE 10260000
+C 10270000
+C### CAUTION ### IF XI,ET,Q ARE SUFFICIENTLY SMALL, THEY ARE SET TO ZER010280000
+C 10290000
+ COMMON /C2/XI2,ET2,Q2,R,R2,R3,R5,Y,D,TT,ALX,ALE,X11,Y11,X32,Y32, 10300000
+ * EY,EZ,FY,FZ,GY,GZ,HY,HZ 10310000
+ DATA F0,F1,F2,EPS/0.D0,1.D0,2.D0,1.D-6/ 10320000
+C----- 10330000
+ IF(DABS(XI).LT.EPS) XI=F0 10340000
+ IF(DABS(ET).LT.EPS) ET=F0 10350000
+ IF(DABS( Q).LT.EPS) Q=F0 10360000
+ XI2=XI*XI 10370000
+ ET2=ET*ET 10380000
+ Q2=Q*Q 10390000
+ R2=XI2+ET2+Q2 10400000
+ R =DSQRT(R2) 10410000
+ IF(R.EQ.F0) RETURN 10420000
+ R3=R *R2 10430000
+ R5=R3*R2 10440000
+ Y =ET*CD+Q*SD 10450000
+ D =ET*SD-Q*CD 10460000
+C----- 10470000
+ IF(Q.EQ.F0) THEN 10480000
+ TT=F0 10490000
+ ELSE 10500000
+ TT=DATAN(XI*ET/(Q*R)) 10510000
+ ENDIF 10520000
+C----- 10530000
+ IF(XI.LT.F0 .AND. Q.EQ.F0 .AND. ET.EQ.F0) THEN 10540002
+ ALX=-DLOG(R-XI) 10550000
+ X11=F0 10560000
+ X32=F0 10570000
+ ELSE 10580000
+ RXI=R+XI 10590002
+ ALX=DLOG(RXI) 10600000
+ X11=F1/(R*RXI) 10610000
+ X32=(R+RXI)*X11*X11/R 10620002
+ ENDIF 10630000
+C----- 10640000
+ IF(ET.LT.F0 .AND. Q.EQ.F0 .AND. XI.EQ.F0) THEN 10650002
+ ALE=-DLOG(R-ET) 10660000
+ Y11=F0 10670000
+ Y32=F0 10680000
+ ELSE 10690000
+ RET=R+ET 10700002
+ ALE=DLOG(RET) 10710000
+ Y11=F1/(R*RET) 10720000
+ Y32=(R+RET)*Y11*Y11/R 10730002
+ ENDIF 10740000
+C----- 10750000
+ EY=SD/R-Y*Q/R3 10760000
+ EZ=CD/R+D*Q/R3 10770000
+ FY=D/R3+XI2*Y32*SD 10780000
+ FZ=Y/R3+XI2*Y32*CD 10790000
+ GY=F2*X11*SD-Y*Q*X32 10800000
+ GZ=F2*X11*CD+D*Q*X32 10810000
+ HY=D*Q*X32+XI*Q*Y32*SD 10820000
+ HZ=Y*Q*X32+XI*Q*Y32*CD 10830000
+ RETURN 10840000
+ END 10850000
Added: cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/dc3d.h
===================================================================
--- cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/dc3d.h 2007-06-26 16:40:29 UTC (rev 7497)
+++ cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/dc3d.h 2007-06-26 16:41:33 UTC (rev 7498)
@@ -0,0 +1,16 @@
+#ifndef __DC3D_H__
+#define __DC3D_H__
+
+void dc3d_(double *alpha,
+ double *x, double *y, double *z,
+ double *depth, double *dip,
+ double *al1, double *al2,
+ double *aw1, double *aw2,
+ double *disl1, double *disl2, double *disl3,
+ double *ux, double *uy, double *uz,
+ double *uxx, double *uyx, double *uzx,
+ double *uxy, double *uyy, double *uzy,
+ double *uxz, double *uyz, double *uzz,
+ double *iret);
+
+#endif
Added: cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/disloc3d.c
===================================================================
--- cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/disloc3d.c 2007-06-26 16:40:29 UTC (rev 7497)
+++ cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/disloc3d.c 2007-06-26 16:41:33 UTC (rev 7498)
@@ -0,0 +1,191 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "dc3d.h"
+#include "disloc3d.h"
+
+//
+// TODO: change flagout,flagout2 to int arrays
+// TODO: flagout2 currently disabled
+//
+
+void disloc3d(double *models, int nmod,
+ double *stations, int nstat,
+ double mu, double nu,
+ double *uout, double *dout, double *sout,
+ double *flagout, double *flagout2)
+{
+ // input arguments
+ // models [nmod x 10]
+ // stations [nstat x 3]
+ // mu [1 x 1]
+ // nu [1 x 1]
+
+ // matrices for return arguments
+ // uout [nstat x 3]
+ // dout [nstat x 9]
+ // sout [nstat x 6]
+ // flagout [nstat x 1]
+ // flagout2 [nstat x nmod]
+
+ double lambda;
+ double theta;
+
+ double *model;
+ double *stat;
+ double *u, *d, *s;
+ double *flag;
+ //double *flag2;
+ double iret;
+
+ double ux, uy, uz;
+ double uxx, uxy, uxz;
+ double uyx, uyy, uyz;
+ double uzx, uzy, uzz;
+
+ double uxt, uyt, uzt;
+ double uxxt, uxyt, uxzt;
+ double uyxt, uyyt, uyzt;
+ double uzxt, uzyt, uzzt;
+
+ double alpha;
+ double x, y, z;
+
+ double strike;
+ double cs,ss;
+
+ double dip;
+ double cd,sd;
+
+ double depth;
+ double disl1, disl2, disl3;
+ double al1, al2;
+ double aw1, aw2;
+
+ int i, j;
+
+
+ lambda = 2*mu*nu/(1 - 2*nu);
+ alpha = (lambda + mu)/(lambda + 2*mu);
+
+ for (i = 0; i < nstat; i++)
+ {
+ stat = &stations[3*i];
+ flag = &flagout[i];
+ //flag2 = &flagout2[nmod*i];
+
+ *flag = 0;
+
+ if (stat[2] > 0)
+ {
+ *flag += 1;
+ fprintf(stderr, "Warning: Positive depth given. (station %d)\n", i);
+ }
+ else
+ {
+ uxt = uyt = uzt = 0;
+ uxxt = uxyt = uxzt = 0;
+ uyxt = uyyt = uyzt = 0;
+ uzxt = uzyt = uzzt = 0;
+
+ for (j = 0; j < nmod; j++)
+ {
+ model = &models[10*j];
+
+ strike = (model[4] - 90)*(M_PI/180.0);
+ cs = cos(strike);
+ ss = sin(strike);
+ dip = model[3];
+ cd = cos(dip * M_PI/180.0);
+ sd = sin(dip * M_PI/180.0);
+ disl1 = model[7];
+ disl2 = model[8];
+ disl3 = model[9];
+ depth = model[2] - 0.5*model[1]*sd;
+ al1 = model[0]/2;
+ al2 = al1;
+ aw1 = model[1]/2;
+ aw2 = aw1;
+ x = cs * (-model[5] + stat[0]) - ss * (-model[6] + stat[1]);
+ y = ss * (-model[5] + stat[0]) + cs * (-model[6] + stat[1]) - 0.5 * cd * model[1];
+ z = stat[2];
+
+ // generate warnings for unphysical models
+ if ((model[2] - sd*model[1] < 0) ||
+ (model[0] <= 0) ||
+ (model[1] <= 0) ||
+ (model[2] < 0))
+ {
+ *flag += 1;
+ fprintf(stderr, "Warning: Unphysical model (station %d, subfault %d)\n",i,j);
+ }
+ else
+ {
+ dc3d_(&alpha, &x, &y, &z, &depth, &dip,
+ &al1, &al2, &aw1, &aw2, &disl1, &disl2, &disl3,
+ &ux, &uy, &uz,
+ &uxx, &uyx, &uzx,
+ &uxy, &uyy, &uzy,
+ &uxz, &uyz, &uzz,
+ &iret);
+
+ // fill flags - the rows contain flags for each station, one row per model/fault
+ *flag += iret;
+ //flag2[j] = iret;
+
+ //
+ // rotate then add
+ //
+ uxt += cs*ux + ss*uy;
+ uyt += -ss*ux + cs*uy;
+ uzt += uz;
+
+ uxxt += cs*cs*uxx + cs*ss*(uxy + uyx) + ss*ss*uyy; // aaa
+ uxyt += cs*cs*uxy - ss*ss*uyx + cs*ss*(-uxx + uyy); // bbb
+ uxzt += cs*uxz + ss*uyz; // ccc
+
+ uyxt += -(ss*(cs*uxx + ss*uxy)) + cs*(cs*uyx + ss*uyy); // ddd
+ uyyt += ss*ss*uxx - cs*ss*(uxy + uyx) + cs*cs*uyy; // eee
+ uyzt += -(ss*uxz) + cs*uyz; // fff
+
+ uzxt += cs*uzx + ss*uzy; // ggg
+ uzyt += -(ss*uzx) + cs*uzy; // hhh
+ uzzt += uzz; // iii
+ }
+ }
+
+ // jump to appropriate offset
+ u = &uout[3*i];
+ d = &dout[9*i];
+ s = &sout[6*i];
+
+ // assign outputs
+ u[0] = uxt;
+ u[1] = uyt;
+ u[2] = uzt;
+
+ d[0] = uxxt;
+ d[1] = uxyt;
+ d[2] = uxzt;
+
+ d[3] = uyxt;
+ d[4] = uyyt;
+ d[5] = uyzt;
+
+ d[6] = uzxt;
+ d[7] = uzyt;
+ d[8] = uzzt;
+
+ // calculate stresses
+ theta = d[0] + d[4] + d[8];
+ s[0] = lambda*theta + 2*mu*d[0];
+ s[1] = mu*(d[1] + d[3]);
+ s[2] = mu*(d[2] + d[6]);
+ s[3] = lambda*theta + 2*mu*d[4];
+ s[4] = mu*(d[5] + d[7]);
+ s[5] = lambda*theta + 2*mu*d[8];
+
+ } // loop over models (subfaults)
+
+ } // loop over stations
+}
Added: cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/disloc3d.h
===================================================================
--- cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/disloc3d.h 2007-06-26 16:40:29 UTC (rev 7497)
+++ cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/disloc3d.h 2007-06-26 16:41:33 UTC (rev 7498)
@@ -0,0 +1,10 @@
+#ifndef __DISLOC3D_H__
+#define __DISLOC3D_H__
+
+void disloc3d(double *models, int nmod,
+ double *stations, int nstat,
+ double mu, double nu,
+ double *uout, double *dout, double *sout,
+ double *flag, double *flag2);
+
+#endif
Added: cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/h5array.c
===================================================================
--- cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/h5array.c 2007-06-26 16:40:29 UTC (rev 7497)
+++ cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/h5array.c 2007-06-26 16:41:33 UTC (rev 7498)
@@ -0,0 +1,95 @@
+#include <stdlib.h>
+#include "hdf5.h"
+#include "h5array.h"
+
+int h5array_read(hid_t loc_id, const char *name,
+ double **data, int *m, int *n)
+{
+ hid_t dataset_id;
+ hid_t dataspace_id;
+ hsize_t dims[2];
+ herr_t status;
+
+ int rank;
+
+
+ dataset_id = H5Dopen(loc_id, name);
+ if (dataset_id < 0)
+ {
+ return -1;
+ }
+
+ dataspace_id = H5Dget_space(dataset_id);
+ rank = H5Sget_simple_extent_ndims(dataspace_id);
+ if (rank != 2)
+ {
+ H5Dclose(dataset_id);
+ H5Sclose(dataspace_id);
+ return -2;
+ }
+
+ status = H5Sget_simple_extent_dims(dataspace_id, dims, NULL);
+ status = H5Sclose(dataspace_id);
+
+ *m = dims[0];
+ *n = dims[1];
+
+ *data = (double *)malloc((dims[0] * dims[1]) * sizeof(double));
+ if (*data == NULL)
+ {
+ H5Dclose(dataset_id);
+ return -3;
+ }
+
+ status = H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
+ H5P_DEFAULT, *data);
+ if (status < 0)
+ {
+ H5Dclose(dataset_id);
+ return -4;
+ }
+
+ status = H5Dclose(dataset_id);
+
+ return 0;
+}
+
+int h5array_write(hid_t loc_id, const char *name,
+ double *data, int m, int n)
+{
+ hsize_t dims[2];
+ hid_t dataspace_id;
+ hid_t dataset_id;
+ herr_t status;
+
+ dims[0] = m;
+ dims[1] = n;
+ dataspace_id = H5Screate_simple(2, dims, NULL);
+ if (dataspace_id < 0)
+ {
+ return -1;
+ }
+
+ dataset_id = H5Dcreate(loc_id, name, H5T_NATIVE_DOUBLE, dataspace_id,
+ H5P_DEFAULT);
+ if (dataset_id < 0)
+ {
+ H5Sclose(dataspace_id);
+ return -2;
+ }
+
+ status = H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
+ H5P_DEFAULT, data);
+ if (status < 0)
+ {
+ H5Sclose(dataspace_id);
+ H5Dclose(dataset_id);
+ return -3;
+ }
+
+ H5Sclose(dataspace_id);
+ H5Dclose(dataset_id);
+
+ return 0;
+}
+
Added: cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/h5array.h
===================================================================
--- cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/h5array.h 2007-06-26 16:40:29 UTC (rev 7497)
+++ cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/h5array.h 2007-06-26 16:41:33 UTC (rev 7498)
@@ -0,0 +1,11 @@
+#ifndef __H5ARRAY_H__
+#define __H5ARRAY_H__
+#include "hdf5.h"
+
+int h5array_read(hid_t loc_id, const char *name,
+ double **data, int *m, int *n);
+
+int h5array_write(hid_t loc_id, const char *name,
+ double *data, int m, int n);
+
+#endif
Added: cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/okadasoln_ng.c
===================================================================
--- cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/okadasoln_ng.c 2007-06-26 16:40:29 UTC (rev 7497)
+++ cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/okadasoln_ng.c 2007-06-26 16:41:33 UTC (rev 7498)
@@ -0,0 +1,502 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+#include "cervelli.h"
+#include "txtarray.h"
+#include "h5array.h"
+
+
+static void linspace(double **x, double a, double b, int n)
+{
+ int i;
+ double dx = (b - a)/n;
+ double *y;
+
+ y = (double *)malloc(n * sizeof(double));
+
+ for (i = 0; i < n; i++)
+ {
+ y[i] = a + i*dx;
+ }
+
+ *x = y;
+}
+
+
+void strike_slip_ng(double *fault_info, int nfaults)
+{
+ int i;
+ double *n;
+ double taperd;
+
+ double *fi;
+ double fx1, fy1;
+ double fx2, fy2;
+ double dip, D, Pr, bd;
+ double ss, ds, ts;
+
+ //
+ // define number of faults to taper with
+ //
+ taperd = 10.0; // (16000.0 - 12000.0)/400
+ linspace(&n, 12000.0, 16000.0, nfaults);
+
+ for (i = 0; i < nfaults; i++)
+ {
+ fi = &fault_info[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/nfaults;
+ 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;
+ }
+
+ free(n);
+}
+
+void reverse_slip_ng(double *fault_info, int nfaults)
+{
+ int i;
+ double *n;
+ double taperd;
+
+ double *fi;
+ double fx1, fy1;
+ double fx2, fy2;
+ double dip, D, Pr, bd;
+ double ss, ds, ts;
+
+ //
+ // define number of faults to taper with
+ //
+ taperd = 10.0; // (16000.0 - 12000.0)/400
+ linspace(&n, 12000.0, 16000.0, nfaults);
+
+ for (i = 0; i < nfaults; i++)
+ {
+ fi = &fault_info[11*i];
+
+ fx1 = 4000.0;
+ fy1 = -n[i];
+
+ fx2 = 4000.0;
+ fy2 = n[i];
+
+ dip = 90.0;
+ D = n[i];
+ Pr = 0.25;
+ bd = 0;
+ ss = 0;
+ ds = 1.0/nfaults;
+ 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;
+ }
+
+ free(n);
+}
+
+
+/*
+ * Note: Be careful with flag2. For a model with
+ * 400 subfaults and 918980 stations, this function would try
+ * to allocate
+ *
+ * (918980*(3+3+9+6+1+400) + 400*(10+11)) * 8 bytes = 2958.8 MiB
+ *
+ * disabling flag2 gives us
+ *
+ * (918980*(3+3+9+6+1) + 400*(10+11)) * 8 bytes = 154.3 MiB
+ *
+ */
+void okada_alloc(int nmod, int nstat,
+ double **fault_info,
+ double **models,
+ double **stations,
+ double **displacement,
+ double **derivatives,
+ double **stress,
+ double **flag,
+ double **flag2)
+{
+ if (fault_info != NULL)
+ *fault_info = (double *)malloc((nmod * 11) * sizeof(double));
+
+ if (models != NULL)
+ *models = (double *)malloc((nmod * 10) * sizeof(double));
+
+ if (stations != NULL)
+ *stations = (double *)malloc((nstat * 3) * sizeof(double));
+
+ if (displacement != NULL)
+ *displacement = (double *)malloc((nstat * 3) * sizeof(double));
+
+ if (derivatives != NULL)
+ *derivatives = (double *)malloc((nstat * 9) * sizeof(double));
+
+ if (stress != NULL)
+ *stress = (double *)malloc((nstat * 6) * sizeof(double));
+
+ if (flag != NULL)
+ *flag = (double *)malloc((nstat * 1) * sizeof(double));
+
+ if (flag2 != NULL)
+ *flag2 = (double *)malloc((nstat * nmod) * sizeof(double));
+}
+
+void okada_free(double *fault_info, double *models, double *stations,
+ double *displacement, double *derivatives, double *stress,
+ double *flag, double *flag2)
+{
+ if (fault_info != NULL) free(fault_info);
+ if (models != NULL) free(models);
+ if (stations != NULL) free(stations);
+ if (displacement != NULL) free(displacement);
+ if (derivatives != NULL) free(derivatives);
+ if (stress != NULL) free(stress);
+ if (flag != NULL) free(flag);
+ if (flag2 != NULL) free(flag2);
+}
+
+
+
+int txt_main(char *infile, char *outfile)
+{
+ int nsubfaults;
+ double *subfaults;
+ double *models;
+
+ int nstations;
+ double *stations;
+ double *displacement;
+ double *derivatives;
+ double *stress;
+
+ double *flag;
+ //double *flag2;
+
+
+ // elastic parameters;
+ double mu = 0.25;
+ double nu = 0.25;
+
+
+
+ /* determine number of subfaults */
+ nsubfaults = 400;
+
+
+ /* read station coordinates */
+ txtarray_read_stations(infile, &stations, &nstations);
+
+ printf("Found %d stations\n", nstations);
+
+ /*
+ * Allocate memory. Skip stations array, since it has
+ * already been allocated above. Also, do not allocate
+ * flag2 array, since it is typically too large.
+ */
+ okada_alloc(nsubfaults, nstations,
+ &subfaults, &models,
+ NULL, &displacement, &derivatives, &stress,
+ &flag, NULL);
+
+
+ /*
+ * Pick one of the following models
+ */
+ strike_slip_ng(subfaults, nsubfaults);
+ //reverse_slip_ng(subfaults, nsubfaults);
+
+
+ /*
+ * Perform the actual calculation
+ */
+ calc_disp_cervelli(mu, nu,
+ models, subfaults, nsubfaults,
+ stations, nstations,
+ displacement, derivatives, stress,
+ flag, NULL);
+
+
+ /*
+ * Write all data to a text file
+ */
+ txtarray_write_all(outfile,
+ "sx sy sz "
+ "ux uy uz "
+ "uxx uxy uxz uyx uyy uyz uzx uzy uzz "
+ "sxx syy szz sxy sxz syz",
+ nstations,
+ stations,
+ displacement,
+ derivatives,
+ stress);
+
+
+ /* free the allocated memory */
+ okada_free(subfaults, models, stations,
+ displacement, derivatives, stress,
+ flag, NULL);
+
+ return 0;
+}
+
+
+
+int split_filepath(const char *filepath, char *filename, char *path, int maxlen)
+{
+ char *x, *y;
+ int i;
+
+ y = strchr(filepath, ':');
+
+ if (y == NULL)
+ {
+ strncpy(filename, filepath, maxlen);
+ path[0] = '\0';
+ return 1;
+ }
+
+ for (i = 0, x = (char *)filepath; (x != y) && (i < maxlen); i++, x++)
+ {
+ filename[i] = *x;
+ }
+ filename[i] = '\0';
+
+ for(i = 0, x = y+1; (*x != '\0') && (i < maxlen); i++, x++)
+ {
+ path[i] = *x;
+ }
+ path[i] = '\0';
+
+ return 0;
+}
+
+int h5_main(char *filepath1, char *filepath2)
+{
+ hid_t file_id;
+ //hid_t group_id;
+ herr_t status;
+
+ char infile[150];
+ char outfile[150];
+ char stations_path[150];
+ char output_group[150];
+
+ int nsubfaults;
+ double *subfaults;
+ double *models;
+
+ int nstations;
+ double *stations;
+ double *displacement;
+ double *derivatives;
+ double *stress;
+
+ double *flag;
+ //double flag2;
+
+ int components;
+
+ // elastic parameters
+ double mu = 0.25;
+ double nu = 0.25;
+
+
+ /* determine number of subfaults */
+ nsubfaults = 400;
+
+
+ /*
+ * Process filepaths
+ */
+
+ split_filepath(filepath1, infile, stations_path, 150);
+ if ('\0' == stations_path[0])
+ {
+ fprintf(stderr, "Error: Could not split file and path from \"%s\"\n", filepath1);
+ return -1;
+ }
+
+ split_filepath(filepath2, outfile, output_group, 150);
+ if (output_group[0] != '\0')
+ {
+ fprintf(stderr, "Error: Writing to an output group currently not supported\n");
+ return -2;
+ }
+
+
+ /*
+ * Read station coordinates
+ */
+
+
+ file_id = H5Fopen(infile, H5F_ACC_RDONLY, H5P_DEFAULT);
+ if (file_id < 0)
+ {
+ fprintf(stderr, "Error: Could not open file %s\n", infile);
+ return -3;
+ }
+
+ status = h5array_read(file_id, stations_path, &stations, &nstations, &components);
+ if (status < 0)
+ {
+ H5Fclose(file_id);
+ return -4;
+ }
+ if (components != 3)
+ {
+ return -5;
+ }
+
+ status = H5Fclose(file_id);
+
+
+ /*
+ * Allocate memory. Skip stations array, since it has
+ * already been allocated above. Also, do not allocate
+ * flag2 array, since it is typically too large.
+ */
+ okada_alloc(nsubfaults, nstations,
+ &subfaults, &models,
+ NULL, &displacement, &derivatives, &stress,
+ &flag, NULL);
+
+ /*
+ * Pick one of the following models
+ */
+ strike_slip_ng(subfaults, nsubfaults);
+ //reverse_slip_ng(subfaults, nsubfaults);
+
+
+ /*
+ * Perform the actual calculation
+ */
+ calc_disp_cervelli(mu, nu,
+ models, subfaults, nsubfaults,
+ stations, nstations,
+ displacement, derivatives, stress,
+ flag, NULL);
+
+
+ /*
+ * Write all data to an HDF5 file (to root group for now)
+ *
+ * TODO: Specify the parent group of these datasets.
+ * Essentially, I need a function that creates
+ * all the necessary groups in a path.
+ */
+
+ file_id = H5Fcreate(outfile, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
+ if (file_id < 0)
+ {
+ H5Fclose(file_id);
+ fprintf(stderr, "Error: Could not write to file %s\n", outfile);
+ return -6;
+ }
+
+ status = h5array_write(file_id, "/stations", stations, nstations, 3);
+ if (status < 0)
+ {
+ H5Fclose(file_id);
+ fprintf(stderr, "Error: Could not write dataset %s\n", "/stations");
+ }
+
+ status = h5array_write(file_id, "/displacement", displacement, nstations, 3);
+ if (status < 0)
+ {
+ H5Fclose(file_id);
+ fprintf(stderr, "Error: Could not write dataset %s\n", "/displacement");
+ }
+
+ status = h5array_write(file_id, "/derivatives", derivatives, nstations, 9);
+ if (status < 0)
+ {
+ H5Fclose(file_id);
+ fprintf(stderr, "Error: Could not write dataset %s\n", "/derivatives");
+ }
+
+ status = h5array_write(file_id, "/stress", stress, nstations, 6);
+ if (status < 0)
+ {
+ H5Fclose(file_id);
+ fprintf(stderr, "Error: Could not write dataset %s\n", "/stress");
+ }
+
+ status = H5Fclose(file_id);
+
+
+ /* free the allocated memory */
+ okada_free(subfaults, models, stations,
+ displacement, derivatives, stress,
+ flag, NULL);
+
+ return 0;
+}
+
+
+/*
+int main(int argc, char *argv[])
+{
+ char arg1[300];
+ char arg2[300];
+
+ if (argc < 4)
+ {
+ fprintf(stderr, "Usage: %s [hdf5 | txt] points.in points.out\n", argv[0]);
+ return -1;
+ }
+
+ snprintf(arg1, 300, argv[2]); // points.in or points.h5:/okada/stations
+ snprintf(arg2, 300, argv[3]); // points.out or data.h5
+
+ if (strcmp(argv[1], "hdf5") == 0)
+ {
+ return h5_main(arg1, arg2);
+ }
+ else if (strcmp(argv[1], "txt") == 0)
+ {
+ return txt_main(arg1, arg2);
+ }
+ else
+ {
+ fprintf(stderr, "Error: Could not understand option \"%s\"\n", argv[3]);
+ }
+
+
+ return EXIT_FAILURE;
+}
+*/
Added: cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/txtarray.c
===================================================================
--- cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/txtarray.c 2007-06-26 16:40:29 UTC (rev 7497)
+++ cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/txtarray.c 2007-06-26 16:41:33 UTC (rev 7498)
@@ -0,0 +1,257 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include "txtarray.h"
+
+
+int txtarray_read1(const char *filename,
+ double **data, int *m, int *n)
+{
+ FILE *fp;
+ char linebuf[500]; // 500 > 14 * (3 + 3 + 9 + 6)
+
+ int i,j;
+ double *x;
+
+ fp = fopen(filename, "r");
+ if (fp == NULL)
+ {
+ return -1;
+ }
+
+ if (fscanf(fp, "# %d %d\n", m, n) != 2)
+ {
+ return -2;
+ }
+
+ fgets(linebuf, 500, fp); // read header line
+
+ *data = (double *)malloc(((*m)*(*n))*sizeof(double));
+
+ for (i = 0; i < *m; i++)
+ {
+ x = &(*data)[(*n) * i];
+ for (j = 0; j < *n; j++)
+ {
+ fscanf(fp, "%lf", &x[j]);
+ }
+ }
+
+ fclose(fp);
+
+ return 0;
+}
+
+
+int txtarray_read_stations(const char *filename,
+ double **stations,
+ int *nstations)
+{
+ int three;
+ int ret;
+
+ ret = txtarray_read1(filename, stations, nstations, &three);
+
+ if (three != 3)
+ {
+ return -1;
+ }
+
+ return ret;
+}
+
+
+/*
+ * TODO: can this family of functions be written using varargs?
+ *
+ * txtarray_write(const char *filename,
+ * const char *header,
+ * int m,
+ * double *data1, int n1,
+ * double *data2, int n2,
+ * ...);
+ *
+ */
+
+int txtarray_write1(const char *filename, const char *header,
+ double *data, int m, int n)
+{
+ FILE *fp;
+ int i,j;
+ double *x;
+
+ fp = fopen(filename, "w");
+ if (fp == NULL)
+ {
+ return -1;
+ }
+
+ fprintf(fp, "# %d %d\n", m, n);
+
+ if (header != NULL)
+ fprintf(fp, "# %s\n", header);
+ else
+ fprintf(fp, "# header\n");
+
+ for (i = 0; i < m; i++)
+ {
+ x = &data[n*i];
+ for (j = 0; j < n; j++)
+ {
+ fprintf(fp, " %0.12e", x[j]);
+ }
+ fprintf(fp, "\n");
+ }
+ fclose(fp);
+
+ return 0;
+}
+
+int txtarray_write2(const char *filename,
+ const char *header,
+ int m,
+ double *data1, int n1,
+ double *data2, int n2)
+{
+ FILE *fp;
+ int i,j;
+ double *x;
+ double *y;
+
+ fp = fopen(filename, "w");
+ if (fp == NULL)
+ {
+ return -1;
+ }
+
+ fprintf(fp, "# %d %d %d\n", m, n1, n2);
+
+ if (header != NULL)
+ fprintf(fp, "# %s\n", header);
+ else
+ fprintf(fp, "# header\n");
+
+ for (i = 0; i < m; i++)
+ {
+ x = &data1[n1 * i];
+ for (j = 0; j < n1; j++)
+ {
+ fprintf(fp, " %0.12e", x[j]);
+ }
+
+ y = &data2[n2 * i];
+ for (j = 0; j < n2; j++)
+ {
+ fprintf(fp, " %0.12e", y[j]);
+ }
+
+ fprintf(fp, "\n");
+ }
+ fclose(fp);
+
+ return 0;
+}
+
+int txtarray_write4(const char *filename,
+ const char *header,
+ int m,
+ double *data1, int n1,
+ double *data2, int n2,
+ double *data3, int n3,
+ double *data4, int n4)
+{
+ FILE *fp;
+ int i,j;
+ double *x, *y, *z, *w;
+
+ fp = fopen(filename, "w");
+ if (fp == NULL)
+ {
+ return -1;
+ }
+
+ fprintf(fp, "# %d %d %d %d %d\n", m, n1, n2, n3, n4);
+
+ if (header != NULL)
+ fprintf(fp, "# %s\n", header);
+ else
+ fprintf(fp, "# header\n");
+
+ for (i = 0; i < m; i++)
+ {
+ x = &data1[n1 * i];
+ for (j = 0; j < n1; j++)
+ {
+ fprintf(fp, " %0.12e", x[j]);
+ }
+
+ y = &data2[n2 * i];
+ for (j = 0; j < n2; j++)
+ {
+ fprintf(fp, " %0.12e", y[j]);
+ }
+
+ z = &data3[n3 * i];
+ for (j = 0; j < n3; j++)
+ {
+ fprintf(fp, " %0.12e", z[j]);
+ }
+
+ w = &data4[n4 * i];
+ for (j = 0; j < n4; j++)
+ {
+ fprintf(fp, " %0.12e", w[j]);
+ }
+
+ fprintf(fp, "\n");
+ }
+ fclose(fp);
+
+ return 0;
+}
+
+int txtarray_write_disp(const char *filename,
+ const char *header,
+ int nstations,
+ double *stations,
+ double *displacements)
+{
+ return txtarray_write2(filename, header, nstations,
+ stations, 3, displacements, 3);
+}
+
+int txtarray_write_deriv(const char *filename,
+ const char *header,
+ int nstations,
+ double *stations,
+ double *derivatives)
+{
+ return txtarray_write2(filename, header, nstations,
+ stations, 3, derivatives, 9);
+}
+
+int txtarray_write_stress(const char *filename,
+ const char *header,
+ int nstations,
+ double *stations,
+ double *stress)
+{
+ return txtarray_write2(filename, header, nstations,
+ stations, 3, stress, 6);
+}
+
+int txtarray_write_all(const char *filename,
+ const char *header,
+ int nstations,
+ double *stations,
+ double *displacements,
+ double *derivatives,
+ double *stress)
+{
+ return txtarray_write4(filename, header,
+ nstations,
+ stations, 3,
+ displacements, 3,
+ derivatives, 9,
+ stress, 6);
+}
+
Added: cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/txtarray.h
===================================================================
--- cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/txtarray.h 2007-06-26 16:40:29 UTC (rev 7497)
+++ cs/cigma/branches/cigma-0.9/benchmarks/disloc3d/txtarray.h 2007-06-26 16:41:33 UTC (rev 7498)
@@ -0,0 +1,56 @@
+#ifndef __TXTARRAY_H__
+#define __TXTARRAY_H__
+
+
+int txtarray_read1(const char *filename,
+ double **data, int *m, int *n);
+
+int txtarray_read_stations(const char *filename,
+ double **stations,
+ int *nstations);
+
+
+int txtarray_write1(const char *filename, const char *header,
+ double *data, int m, int n);
+
+int txtarray_write2(const char *filename,
+ const char *header,
+ int m,
+ double *data1, int n1,
+ double *data2, int n2);
+
+int txtarray_write4(const char *filename,
+ const char *header,
+ int m,
+ double *data1, int n1,
+ double *data2, int n2,
+ double *data3, int n3,
+ double *data4, int n4);
+
+int txtarray_write_disp(const char *filename,
+ const char *header,
+ int nstations,
+ double *stations,
+ double *displacements);
+
+int txtarray_write_deriv(const char *filename,
+ const char *header,
+ int nstations,
+ double *stations,
+ double *derivatives);
+
+int txtarray_write_stress(const char *filename,
+ const char *header,
+ int nstations,
+ double *stations,
+ double *displacements);
+
+int txtarray_write_all(const char *filename,
+ const char *header,
+ int nstations,
+ double *stations,
+ double *displacements,
+ double *derivatives,
+ double *stress);
+
+#endif
More information about the cig-commits
mailing list