elamodel.c

su 的源代码库

/* Copyright (c) Colorado School of Mines, 2006.*/
/* All rights reserved.                       */

/* Copyright (c) Colorado School of Mines, 1999.*/
/* All rights reserved.                       */

/* ELAMODEL: $Test Release: 1.1 $ ; $Date: 1997/03/05 16:45:08 $	*/

#include "par.h"
#include "tri.h"
#include "elastic.h"

/*********************** self documentation **********************/
char *sdoc[] = {
"									",
" ELAMODEL - make piecewise homogeneous anisotropic model    		",
"									",
" elamodel >modelfile fill= [optional parameters]   			",
" 									",
" Input Parameters: 							",
" xmin=0.0               minimum horizontal coordinate (x) 		",
" xmax=1.0               maximum horizontal coordinate (x) 		",
" zmin=0.0               minimum vertical coordinate (z) 		",
" zmax=1.0               maximum vertical coordinate (z) 		",
" xedge=                 x coordinates of an edge 			",
" zedge=                 z coordinates of an edge 			",
" kedge=                 array of indices used to identify edges 	",
" fill=    iso      	 x,z,v_p,v_s,rho   				",
"          tiso      	 x,z,v_p,v_s,epsilon,delta,gamma,phi,rho	",
"          ani           x,z,a1111,a3333,a1133,a1313,a1113,a3313        ",
"                            a1212,a2323,a1223,rho                      ",
" maxangle=5.0           maximum angle (deg) between adjacent edge segments ",
"   									",
" Notes: 								",
" More than set of xedge and zedge parameters may be 		        ",
" specified, but the numbers of these parameters must be equal. 	",
" 									",
" Within each set, vertices will be connected by fixed edges. 		",
" 									",
" Edge indices in the k array are used to identify edges 		",
" specified by the x and z parameters.  The first k index 		",
" corresponds to the first set of x and z parameters, the 		",
" second k index corresponds to the second set, and so on. 		",
" 									",
" After all vertices have been inserted into the model,	the fill        ",
" parameters is used to fill closed regions bounded by fixed edges.     ",
" Three input modes are available:                                      ",
" Isotropic blocks:     x,z,v_p,v_s,rho                                 ",
" Transversely iso:     x,z,v_p,v_s,epsilon,delta,gamma,phi,rho         ",
" General 2D aniso:     x,z,a1111,a3333,a1133,a1313,a1113,a3313         ",
"                       a1212,a2323,a1223,rho                           ",
" Hereby:  ",
" x,z			   coordinates of one point in a bounded region ",
" v_p,v_s		   P, S-wave velocity along symmetry axis       ",
" epsilon, delta, gammma   Thomsen's parameters              ",
" rho 			   density 			     ",
" phi			   angle of symmetry axes with vertical ",
" aijkl			   density normalized stiffness coefficients ",
"  ",
" Each block can be defined by different input modes. The number of     ",
" input parameters defines the input type. Incorrect number of input    ",
" parameters result in an Error-message					",
" 									",
NULL};

/*
 *
 * AUTHOR:  Dave Hale, Colorado School of Mines, 02/12/91
 * modified : Andreas Rueger, Colorado School of Mines, 01/18/94
 * built anisotropic models
 *
 */
/**************** end self doc ***********************************/

/* prototypes for functions defined and used internally */ 

static void smoothTwoSegments (float maxangle, int i,
	int nd, float ud[], float xd[][4], float yd[][4],
	int *n, float *u[], float *x[], float *y[]);

static void smoothLineSegments (float maxangle,	int nin, float *xin,
	float *yin, int *nout, float **xout, float **yout, float **txout,
	float **tyout, float **cout);

static Model *makemod (float xmin, float zmin, float xmax, float zmax,
	int ne, int *ke, int *ns, float **xs, float **zs,
	float **txs, float **tzs, float **cs);

static Vertex *newVertex (Model *m, Vertex *vlast,
	float x, float z);

static void setEdgeAttributes (Vertex *v1, Vertex *v2, int k);

static void setEdgeUseAttributes (Vertex *v1, Vertex *v2,float tx1,
	 float tz1, float c1, float tx2, float tz2, float c2);

static void fillcoeff (Model *m, int nr, float **sra);

static void setcoeff (Tri *t, float a1111, float a3333, float a1133, float
 	a1313, float a1113, float a3313, float a1212, float a2323,
	float a1223, float rho, int mindex);

static void checkfill (Model *m, int nr, float **sra);

int rottensors (float *aijkl, float angle);

/* the main program */
int main (int argc, char **argv)
{
	int nxe,nze,nke,ne,nr,*nxz,*ns,ie,ir,*ke,num,i;
	float xmin,zmin,xmax,zmax,maxangle,temp,
		**xe,**ze,**xs,**zs,**txs,**tzs,**cs,*sr,**sra;
	float angle;
	Model *m;

	/* hook up getpar to handle the parameters */
	initargs(argc,argv);
	requestdoc(0);
	
	/***************** get optional parameters **********************/

	if (!getparfloat("xmin",&xmin)) xmin = 0.0;
	if (!getparfloat("xmax",&xmax)) xmax = 1.0;
	if (!getparfloat("zmin",&zmin)) zmin = 0.0;
	if (!getparfloat("zmax",&zmax)) zmax = 1.0;
	if (!getparfloat("maxangle",&maxangle)) maxangle = 5.0;
	maxangle *= PI/180.0;

	/***************** input interfaces ******************************/

	nxe = countparname("xedge");
	nze = countparname("zedge");

	if (nxe!=nze)
	 err("\n ERROR number of xedge must equal number of zedge\n");

	ne = nxe;
	if (ne>0) {	
		xe = (float**)malloc(ne*sizeof(float*));
		ze = (float**)malloc(ne*sizeof(float*));
		nxz = alloc1int(ne);
	} else {
		xe = (float**)malloc(sizeof(float*));
		ze = (float**)malloc(sizeof(float*));
		nxz = alloc1int(1);
	}	

	for (ie=0; ie<ne; ++ie) {
		nxe = countnparval(ie+1,"xedge");
		nze = countnparval(ie+1,"zedge");
		if (nxe!=nze) 
		 err("\n ERROR number of xedge must equal number of zedge\n");

		nxz[ie] = nxe;
		xe[ie] = alloc1float(nxz[ie]);
		ze[ie] = alloc1float(nxz[ie]);
		getnparfloat(ie+1,"xedge",xe[ie]);
		getnparfloat(ie+1,"zedge",ze[ie]);
	}


	nke = countparval("kedge");
	ke = alloc1int(ne);
	if (nke>ne) err("\n ERROR more <kedge> values than edges specified\n");
	if (nke<ne) err("\n ERROR less <kedge> values than edges specified\n");

	getparint("kedge",ke);


	/*************** read seismic parameters *****************************/
	nr = countparname("fill");
	if(nr ==0)
		 err("\n ERROR required parameters <fill> missing \n");

	sr = alloc1float(12);
	sra = alloc2float(13,nr);

	for (ir=0; ir<nr; ++ir) {

		num=getnparfloat(ir+1,"fill",sr);
		if(num != 5 && num != 9 && num != 12)
		err("\nERROR wrong number of input parameters in <fill> #%i\n",
		 	ir+1);

                /****** isotropic input *******/
		if(num == 5){
			sra[ir][0] = sr[0];
		 	sra[ir][1] = sr[1];
			sra[ir][2] = sra[ir][3] = sr[2]*sr[2];
			sra[ir][5] = sr[3]*sr[3];
 			sra[ir][4] = sra[ir][2] - 2*sra[ir][5];
			sra[ir][6] = sra[ir][7] = 0;
			sra[ir][8] = sra[ir][9] = sra[ir][5];
			sra[ir][10] = 0;
			sra[ir][11] = sr[4];
			sra[ir][12] = 0;

		/*** transversly isotropic ***/
		} else if(num == 9){	
			sra[ir][0] = sr[0];
			sra[ir][1] = sr[1];
			sra[ir][3] = sr[2]*sr[2];
			sra[ir][5] = sr[3]*sr[3];
			sra[ir][2] = sra[ir][3]+2*sr[4]*sra[ir][3];

			temp=2*sr[5]*sra[ir][3]*(sra[ir][3]-sra[ir][5])
			+(sra[ir][3]-sra[ir][5])*(sra[ir][3]-sra[ir][5]);

			if(temp <0) err("ERROR: a1133 not reel \n");

			sra[ir][4] = sqrt(temp)-sra[ir][5];

			sra[ir][6] = sra[ir][7] = sra[ir][10] =0;

			sra[ir][8] = 2*sra[ir][5]*sr[6]  + sra[ir][5];

			sra[ir][9] = sra[ir][5];

			sra[ir][11] = sr[8];
			sra[ir][12] = 1;


      			angle=sr[7]*PI/180.0;

    			if(angle != 0){

 			  /******** do the rotation ***********/
			  if(rottensors(sra[ir], angle) !=1)
				err(" Error \n rotation failed \n");

			}

		} else if (num==12){
			/****** anisotropic ************/

			for(i=0;i<12;i++)
				sra[ir][i] = sr[i];

			sra[ir][12]=2;

		} else {
		  err("\n ERROR wrong input parameters in <fill> #%i\n",ir+1);
		}
	}

	/*************** smooth interfaces *********************************/
	if (ne>0) {
		ns = alloc1int(ne);
 		xs = (float**)malloc(ne*sizeof(float*));
		zs = (float**)malloc(ne*sizeof(float*));
		txs = (float**)malloc(ne*sizeof(float*));
		tzs = (float**)malloc(ne*sizeof(float*));
		cs = (float**)malloc(ne*sizeof(float*));
		for (ie=0; ie<ne; ++ie)
			smoothLineSegments(maxangle,nxz[ie],xe[ie],ze[ie],
			&ns[ie],&xs[ie],&zs[ie],&txs[ie],&tzs[ie],&cs[ie]);
	} else {
		ns = alloc1int(1);
 		xs = (float**)malloc(sizeof(float*));
		zs = (float**)malloc(sizeof(float*));
		txs = (float**)malloc(sizeof(float*));
		tzs = (float**)malloc(sizeof(float*));
		cs = (float**)malloc(sizeof(float*));
	}

	/*************** create the  model *********************************/
	m = makemod(xmin,zmin,xmax,zmax,ne,ke,ns,xs,zs,txs,tzs,cs);

	/********* fill regions with elastic stiffness coefficient *********/
	fillcoeff(m,nr,sra);

	/*************** check if all regions are filled *******************/
	checkfill(m,nr,sra);

	/********************* write model *********************************/
	writeModel(m,stdout);

        return 1;
}

/* build model by connecting vertices of smoothed edges */
static Model *makemod (float xmin, float zmin, float xmax, float zmax,
	int ne, int *ke, int *ns, float **xs, float **zs, 
	float **txs, float **tzs, float **cs)
{
	int ie,is,islast;
	float txs1,tzs1,cs1,txs2,tzs2,cs2;
	Vertex *vlast,*v;
	Model *m;
	
	/* initialize model */
	m = makeModel(zmin,xmin,zmax,xmax);
	m->eps = 1.0e-5*sqrt(pow(xmax-xmin,2.0)+pow(zmax-zmin,2.0));
	m->sfa = sizeof(FaceAttributes);
	m->seua = sizeof(EdgeUseAttributes);
	m->sea = sizeof(EdgeAttributes);

	/* loop over smoothed edges */
	for (ie=0; ie<ne; ++ie) {
		
		/* add first vertex to model */
		v = newVertex(m,NULL,xs[ie][0],zs[ie][0]);
		vlast = v;
		islast = 0;
		
		/* loop over vertices in smoothed edge */
		for (is=1; is<ns[ie]; ++is) {
			
			/* add vertex to model */
			v = newVertex(m,vlast,
				xs[ie][is],zs[ie][is]);
			
			/* if new vertex is too close to last vertex */
			if (v==vlast) continue;
			
			/* connect vertex to last vertex with fixed edge */
			fixEdgeBetweenVertices(vlast,v);
			
			/* set edge attributes */
			setEdgeAttributes(vlast,v,ke[ie]);
			
			/* set edge-use attributes */
			txs1 = txs[ie][islast];
			tzs1 = tzs[ie][islast];
			cs1 = cs[ie][islast];
			txs2 = -txs[ie][is];
			tzs2 = -tzs[ie][is];
			cs2 = -cs[ie][is];
			setEdgeUseAttributes(vlast,v,
				txs1,tzs1,cs1,txs2,tzs2,cs2);
			
			/* remember last vertex */
			vlast = v;
			islast = is;
		}
	}
	
	
	/* return model */
	return m;
}




static void smoothLineSegments (float maxangle,	int nin, float *xin, float *yin,
	int *nout, float **xout, float **yout, float **txout, float **tyout,
	float **cout)
/* smooth line segments and compute tangent vectors and curvatures */
{