[cig-commits] r8963 - in mc/3D/CitcomS/trunk: lib tests
tan2 at geodynamics.org
tan2 at geodynamics.org
Thu Dec 20 12:30:14 PST 2007
Author: tan2
Date: 2007-12-20 12:30:13 -0800 (Thu, 20 Dec 2007)
New Revision: 8963
Added:
mc/3D/CitcomS/trunk/tests/gnomonic.c
Modified:
mc/3D/CitcomS/trunk/lib/Full_sphere_related.c
mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c
mc/3D/CitcomS/trunk/lib/Full_version_dependent.c
mc/3D/CitcomS/trunk/lib/Sphere_harmonics.c
mc/3D/CitcomS/trunk/lib/Sphere_util.c
mc/3D/CitcomS/trunk/lib/global_defs.h
Log:
Generated global mesh by great circles. Tracer module can split the cap in the map view now.
Modified: mc/3D/CitcomS/trunk/lib/Full_sphere_related.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_sphere_related.c 2007-12-20 20:27:41 UTC (rev 8962)
+++ mc/3D/CitcomS/trunk/lib/Full_sphere_related.c 2007-12-20 20:30:13 UTC (rev 8963)
@@ -25,6 +25,8 @@
*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/
+
+
/* Functions relating to the building and use of mesh locations ... */
@@ -33,41 +35,181 @@
#include "element_definitions.h"
#include "global_defs.h"
-void full_coord_of_cap(E,m,icap)
- struct All_variables *E;
- int icap,m;
- {
- int i,j,k,lev,temp,elx,ely,nox,noy,noz,node,nodes;
- int lelx,lely,lnox,lnoy,lnoz;
- double x[5],y[5],z[5],xx[5],yy[5],zz[5];
- double *theta1,*fi1,*theta2,*fi2,*theta,*fi,*SX[2];
+/**************************************************************/
+/* This function transforms theta and phi to new coords */
+/* u and v using gnomonic projection. */
+/* See http://mathworld.wolfram.com/GnomonicProjection.html */
+
+void spherical_to_uv2(double center[2], int len,
+ double *theta, double *phi,
+ double *u, double *v)
+{
+ double theta_f, phi_f;
+ double cos_tf, sin_tf;
+ double cosc, cost, sint, cosp2, sinp2;
+ int i;
+
+ /* theta_f and phi_f are the reference points at the midpoint of the cap */
+
+ theta_f = center[0];
+ phi_f = center[1];
+
+ cos_tf = cos(theta_f);
+ sin_tf = sin(theta_f);
+
+ for(i=0; i<len; i++) {
+ cost = cos(theta[i]);
+ sint = sin(theta[i]);
+
+ cosp2 = cos(phi[i] - phi_f);
+ sinp2 = sin(phi[i] - phi_f);
+
+ cosc = cos_tf * cost + sin_tf * sint * cosp2;
+ cosc = 1.0 / cosc;
+
+ u[i] = sint * sinp2 * cosc;
+ v[i] = (sin_tf * cost - cos_tf * sint * cosp2) * cosc;
+ }
+ return;
+}
+
+
+/**************************************************************/
+/* This function transforms u and v to spherical coord */
+/* theta and phi using inverse gnomonic projection. */
+/* See http://mathworld.wolfram.com/GnomonicProjection.html */
+
+void uv_to_spherical(double center[2], int len,
+ double *u, double *v,
+ double *theta, double *phi)
+{
+ double theta_f, phi_f, cos_tf, sin_tf;
+ double x, y, r, c;
+ double cosc, sinc;
+ int i;
+
+ /* theta_f and phi_f are the reference points at the midpoint of the cap */
+
+ theta_f = center[0];
+ phi_f = center[1];
+
+ cos_tf = cos(theta_f);
+ sin_tf = sin(theta_f);
+
+ for(i=0; i<len; i++) {
+ x = u[i];
+ y = v[i];
+ r = sqrt(x*x + y*y);
+
+ /* special case: r=0, then (u,v) is the reference point */
+ if(r == 0) {
+ theta[i] = theta_f;
+ phi[i] = phi_f;
+ continue;
+ }
+
+ /* c = atan(r); cosc = cos(c); sinc = sin(c);*/
+ cosc = 1.0 / sqrt(1 + r*r);
+ sinc = sqrt(1 - cosc*cosc);
+
+ theta[i] = acos(cosc * cos_tf +
+ y * sinc * sin_tf / r);
+ phi[i] = phi_f + atan(x * sinc /
+ (r * sin_tf * cosc - y * cos_tf * sinc));
+ }
+ return;
+}
+
+
+/* Find the intersection point of two lines */
+/* The lines are: (x[0], y[0]) to (x[1], y[1]) */
+/* (x[2], y[2]) to (x[3], y[3]) */
+/* If found, the intersection point is stored */
+/* in (px, py) and return 1 */
+/* If not found, return 0 */
+
+static int find_intersection(double *x, double *y,
+ double *px, double *py)
+{
+ double a1, b1, c1;
+ double a2, b2, c2;
+ double denom;
+
+ a1 = y[1] - y[0];
+ b1 = x[0] - x[1];
+ c1 = x[1]*y[0] - x[0]*y[1];
+
+ a2 = y[3] - y[2];
+ b2 = x[2] - x[3];
+ c2 = x[3]*y[2] - x[2]*y[3];
+
+ denom = a1*b2 - a2*b1;
+ if (denom == 0) return 0; /* the lines are parallel! */
+
+ *px = (b1*c2 - b2*c1)/denom;
+ *py = (a2*c1 - a1*c2)/denom;
+ return 1;
+}
+
+
+void full_coord_of_cap(struct All_variables *E, int m, int icap)
+{
+ int i, j, k, lev, temp, elx, ely, nox, noy, noz;
+ int node, snode;
+ int lelx, lely, lnox, lnoy, lnoz;
+ int ok;
+ double x[5], y[5], z[5], referencep[2];
+ double *theta0, *fi0;
+ double *tt1, *tt2, *tt3, *tt4, *ff1, *ff2, *ff3, *ff4;
+ double *u1, *u2, *u3, *u4, *v1, *v2, *v3, *v4;
+ double *px, *py, *qx, *qy;
+ double theta, fi, cost, sint, cosf, sinf;
+ double a, b;
double myatan();
void even_divide_arc12();
temp = max(E->mesh.NOY[E->mesh.levmax],E->mesh.NOX[E->mesh.levmax]);
- theta1 = (double *)malloc((temp+1)*sizeof(double));
- fi1 = (double *)malloc((temp+1)*sizeof(double));
- theta2 = (double *)malloc((temp+1)*sizeof(double));
- fi2 = (double *)malloc((temp+1)*sizeof(double));
- theta = (double *)malloc((temp+1)*sizeof(double));
- fi = (double *)malloc((temp+1)*sizeof(double));
+ theta0 = (double *)malloc((temp+1)*sizeof(double));
+ fi0 = (double *)malloc((temp+1)*sizeof(double));
- temp = E->mesh.NOY[E->mesh.levmax]*E->mesh.NOX[E->mesh.levmax]; /* allocate enough for the entire cap */
+ tt1 = (double *)malloc((temp+1)*sizeof(double));
+ tt2 = (double *)malloc((temp+1)*sizeof(double));
+ tt3 = (double *)malloc((temp+1)*sizeof(double));
+ tt4 = (double *)malloc((temp+1)*sizeof(double));
- SX[0] = (double *)malloc((temp+1)*sizeof(double));
- SX[1] = (double *)malloc((temp+1)*sizeof(double));
+ ff1 = (double *)malloc((temp+1)*sizeof(double));
+ ff2 = (double *)malloc((temp+1)*sizeof(double));
+ ff3 = (double *)malloc((temp+1)*sizeof(double));
+ ff4 = (double *)malloc((temp+1)*sizeof(double));
+ u1 = (double *)malloc((temp+1)*sizeof(double));
+ u2 = (double *)malloc((temp+1)*sizeof(double));
+ u3 = (double *)malloc((temp+1)*sizeof(double));
+ u4 = (double *)malloc((temp+1)*sizeof(double));
+
+ v1 = (double *)malloc((temp+1)*sizeof(double));
+ v2 = (double *)malloc((temp+1)*sizeof(double));
+ v3 = (double *)malloc((temp+1)*sizeof(double));
+ v4 = (double *)malloc((temp+1)*sizeof(double));
+
+ temp = E->mesh.NOY[E->mesh.levmax] * E->mesh.NOX[E->mesh.levmax];
+ px = malloc((temp+1)*sizeof(double));
+ py = malloc((temp+1)*sizeof(double));
+ qx = malloc((temp+1)*sizeof(double));
+ qy = malloc((temp+1)*sizeof(double));
+
+ /* 4 corners of the cap */
for (i=1;i<=4;i++) {
x[i] = sin(E->sphere.cap[icap].theta[i])*cos(E->sphere.cap[icap].fi[i]);
y[i] = sin(E->sphere.cap[icap].theta[i])*sin(E->sphere.cap[icap].fi[i]);
z[i] = cos(E->sphere.cap[icap].theta[i]);
- }
-
- for (lev=E->mesh.levmin;lev<=E->mesh.levmax;lev++) {
+ }
+ for (lev=E->mesh.levmin;lev<=E->mesh.levmax;lev++) {
+
elx = E->lmesh.ELX[lev]*E->parallel.nprocx;
ely = E->lmesh.ELY[lev]*E->parallel.nprocy;
nox = elx+1;
@@ -78,67 +220,176 @@
lnox = lelx+1;
lnoy = lely+1;
lnoz = E->lmesh.NOZ[lev];
- /* evenly divide arc linking 1 and 2, and the arc linking 4 and 3 */
- even_divide_arc12(elx,x[1],y[1],z[1],x[2],y[2],z[2],theta1,fi1);
- even_divide_arc12(elx,x[4],y[4],z[4],x[3],y[3],z[3],theta2,fi2);
+ /* evenly divide arc linking 1 and 2 */
+ even_divide_arc12(elx,x[1],y[1],z[1],x[2],y[2],z[2],theta0,fi0);
- for (j=1;j<=nox;j++) {
+ /* pick up only points within this processor */
+ for (j=0, i=E->lmesh.nxs; j<lnox; j++, i++) {
+ tt1[j] = theta0[i];
+ ff1[j] = fi0[i];
+ }
- /* pick up the two ends */
- xx[1] = sin(theta1[j])*cos(fi1[j]);
- yy[1] = sin(theta1[j])*sin(fi1[j]);
- zz[1] = cos(theta1[j]);
- xx[2] = sin(theta2[j])*cos(fi2[j]);
- yy[2] = sin(theta2[j])*sin(fi2[j]);
- zz[2] = cos(theta2[j]);
+ /* evenly divide arc linking 4 and 3 */
+ even_divide_arc12(elx,x[4],y[4],z[4],x[3],y[3],z[3],theta0,fi0);
+ /* pick up only points within this processor */
+ for (j=0, i=E->lmesh.nxs; j<lnox; j++, i++) {
+ tt2[j] = theta0[i];
+ ff2[j] = fi0[i];
+ }
- even_divide_arc12(ely,xx[1],yy[1],zz[1],xx[2],yy[2],zz[2],theta,fi);
+ /* evenly divide arc linking 1 and 4 */
+ even_divide_arc12(ely,x[1],y[1],z[1],x[4],y[4],z[4],theta0,fi0);
- for (k=1;k<=noy;k++) {
- nodes = j + (k-1)*nox;
- SX[0][nodes] = theta[k];
- SX[1][nodes] = fi[k];
- }
+ /* pick up only points within this processor */
+ for (k=0, i=E->lmesh.nys; k<lnoy; k++, i++) {
+ tt3[k] = theta0[i];
+ ff3[k] = fi0[i];
+ }
+ /* evenly divide arc linking 2 and 3 */
+ even_divide_arc12(ely,x[2],y[2],z[2],x[3],y[3],z[3],theta0,fi0);
- } /* end for j */
+ /* pick up only points within this processor */
+ for (k=0, i=E->lmesh.nys; k<lnoy; k++, i++) {
+ tt4[k] = theta0[i];
+ ff4[k] = fi0[i];
+ }
- /* get the coordinates for the entire cap */
+ /* compute the intersection point of these great circles */
+ /* the point is first found in u-v space and project back */
+ /* to theta-phi space later */
- for (j=1;j<=lnox;j++)
- for (k=1;k<=lnoy;k++) {
- nodes = (j+(E->lmesh.NXS[lev]-1))+(k-1+(E->lmesh.NYS[lev]-1))*nox;
- for (i=1;i<=lnoz;i++) {
- node = i + (j-1)*lnoz + (k-1)*lnox*lnoz;
+ referencep[0] = E->sphere.cap[icap].theta[2];
+ referencep[1] = E->sphere.cap[icap].fi[2];
- /* theta,fi,and r coordinates */
- E->SX[lev][m][1][node] = SX[0][nodes];
- E->SX[lev][m][2][node] = SX[1][nodes];
- E->SX[lev][m][3][node] = E->sphere.R[lev][i];
+ spherical_to_uv2(referencep, lnox, tt1, ff1, u1, v1);
+ spherical_to_uv2(referencep, lnox, tt2, ff2, u2, v2);
+ spherical_to_uv2(referencep, lnoy, tt3, ff3, u3, v3);
+ spherical_to_uv2(referencep, lnoy, tt4, ff4, u4, v4);
- /* x,y,and z oordinates */
- E->X[lev][m][1][node] =
- E->sphere.R[lev][i]*sin(SX[0][nodes])*cos(SX[1][nodes]);
- E->X[lev][m][2][node] =
- E->sphere.R[lev][i]*sin(SX[0][nodes])*sin(SX[1][nodes]);
- E->X[lev][m][3][node] =
- E->sphere.R[lev][i]*cos(SX[0][nodes]);
- }
+ snode = 0;
+ for(k=0; k<lnoy; k++) {
+ x[2] = u3[k];
+ y[2] = v3[k];
+
+ x[3] = u4[k];
+ y[3] = v4[k];
+
+ for(j=0; j<lnox; j++) {
+ x[0] = u1[j];
+ y[0] = v1[j];
+
+ x[1] = u2[j];
+ y[1] = v2[j];
+
+ ok = find_intersection(x, y, &a, &b);
+ if(!ok) {
+ fprintf(stderr, "Error(Full_coord_of_cap): cannot find intersection point!\n");
+ exit(10);
}
- } /* end for lev */
+ px[snode] = a;
+ py[snode] = b;
+ snode++;
+ }
+ }
- free ((void *)SX[0]);
- free ((void *)SX[1]);
- free ((void *)theta );
- free ((void *)theta1);
- free ((void *)theta2);
- free ((void *)fi );
- free ((void *)fi1 );
- free ((void *)fi2 );
+ uv_to_spherical(referencep, snode, px, py, qx, qy);
+ /* replace (qx, qy) by (tt?, ff?) for points on the edge */
+ if(E->parallel.me_loc[2] == 0) {
+ /* x boundary */
+ for(k=0; k<lnox; k++) {
+ i = k;
+ qx[i] = tt1[k];
+ qy[i] = ff1[k];
+ }
+ }
+
+ if(E->parallel.me_loc[2] == E->parallel.nprocy-1) {
+ /* x boundary */
+ for(k=0; k<lnox; k++) {
+ i = k + (lnoy - 1) * lnox;
+ qx[i] = tt2[k];
+ qy[i] = ff2[k];
+ }
+ }
+
+ if(E->parallel.me_loc[1] == 0) {
+ /* y boundary */
+ for(k=0; k<lnoy; k++) {
+ i = k * lnox;
+ qx[i] = tt3[k];
+ qy[i] = ff3[k];
+ }
+ }
+
+ if(E->parallel.me_loc[1] == E->parallel.nprocx-1) {
+ /* y boundary */
+ for(k=0; k<lnoy; k++) {
+ i = (k + 1) * lnox - 1;
+ qx[i] = tt4[k];
+ qy[i] = ff4[k];
+ }
+ }
+
+ /* store the node location in spherical and cartesian coordinates */
+ for (k=0; k<snode; k++) {
+ theta = qx[k];
+ fi = qy[k];
+
+ cost = cos(theta);
+ sint = sin(theta);
+ cosf = cos(fi);
+ sinf = sin(fi);
+
+ for (i=1; i<=lnoz; i++) {
+ node = i + k*lnoz;
+
+ /* theta,fi,and r coordinates */
+ E->SX[lev][m][1][node] = theta;
+ E->SX[lev][m][2][node] = fi;
+ E->SX[lev][m][3][node] = E->sphere.R[lev][i];
+
+ /* x,y,and z oordinates */
+ E->X[lev][m][1][node] = E->sphere.R[lev][i]*sint*cosf;
+ E->X[lev][m][2][node] = E->sphere.R[lev][i]*sint*sinf;
+ E->X[lev][m][3][node] = E->sphere.R[lev][i]*cost;
+ }
+ }
+
+ } /* end for lev */
+
+ free ((void *)theta0);
+ free ((void *)fi0 );
+
+ free ((void *)tt1 );
+ free ((void *)tt2 );
+ free ((void *)tt3 );
+ free ((void *)tt4 );
+
+ free ((void *)ff1 );
+ free ((void *)ff2 );
+ free ((void *)ff3 );
+ free ((void *)ff4 );
+
+ free ((void *)u1 );
+ free ((void *)u2 );
+ free ((void *)u3 );
+ free ((void *)u4 );
+
+ free ((void *)v1 );
+ free ((void *)v2 );
+ free ((void *)v3 );
+ free ((void *)v4 );
+
+ free ((void *)px );
+ free ((void *)py );
+ free ((void *)qx );
+ free ((void *)qy );
+
return;
- }
+}
Modified: mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c 2007-12-20 20:27:41 UTC (rev 8962)
+++ mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c 2007-12-20 20:30:13 UTC (rev 8963)
@@ -1834,13 +1834,6 @@
fprintf(E->trace.fpt,"\nTracing Activated! (proc: %d)\n",E->parallel.me);
fprintf(E->trace.fpt," Allen K. McNamara 12-2003\n\n");
-
- if ((E->parallel.nprocx > 1) || (E->parallel.nprocy > 1)) {
- fprintf(E->trace.fpt,"Sorry - Tracer code does not (yet) work if nprocx or nprocy is greater than 1\n");
- fflush(E->trace.fpt);
- parallel_process_termination();
- }
-
if (E->trace.ic_method==0)
{
fprintf(E->trace.fpt,"Generating New Tracer Array\n");
Modified: mc/3D/CitcomS/trunk/lib/Full_version_dependent.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_version_dependent.c 2007-12-20 20:27:41 UTC (rev 8962)
+++ mc/3D/CitcomS/trunk/lib/Full_version_dependent.c 2007-12-20 20:30:13 UTC (rev 8963)
@@ -126,7 +126,40 @@
}
+
/* =================================================
+ rotate the mesh by a rotation matrix
+ =================================================*/
+static void full_rotate_mesh(struct All_variables *E, double dircos[4][4],
+ int m, int icap)
+{
+ int i,lev;
+ double t[4], myatan();
+
+ for (lev=E->mesh.levmin;lev<=E->mesh.levmax;lev++) {
+ for (i=1;i<=E->lmesh.NNO[lev];i++) {
+ t[0] = E->X[lev][m][1][i]*dircos[1][1]+
+ E->X[lev][m][2][i]*dircos[1][2]+
+ E->X[lev][m][3][i]*dircos[1][3];
+ t[1] = E->X[lev][m][1][i]*dircos[2][1]+
+ E->X[lev][m][2][i]*dircos[2][2]+
+ E->X[lev][m][3][i]*dircos[2][3];
+ t[2] = E->X[lev][m][1][i]*dircos[3][1]+
+ E->X[lev][m][2][i]*dircos[3][2]+
+ E->X[lev][m][3][i]*dircos[3][3];
+
+ E->X[lev][m][1][i] = t[0];
+ E->X[lev][m][2][i] = t[1];
+ E->X[lev][m][3][i] = t[2];
+ E->SX[lev][m][1][i] = acos(t[2]/E->SX[lev][m][3][i]);
+ E->SX[lev][m][2][i] = myatan(t[1],t[0]);
+ }
+ } /* lev */
+
+ return;
+}
+
+/* =================================================
Standard node positions including mesh refinement
================================================= */
@@ -136,6 +169,7 @@
{
int i,j,k,ii,lev;
double ro,dr,*rr,*RR,fo;
+ double dircos[4][4];
float tt1;
int step,nn;
char output_file[255], a[255];
@@ -209,30 +243,54 @@
free ((void *) rr);
free ((void *) RR);
- ro = -0.5*(M_PI/4.0)/E->mesh.elx;
- fo = 0.0;
-
- E->sphere.dircos[1][1] = cos(ro)*cos(fo);
- E->sphere.dircos[1][2] = cos(ro)*sin(fo);
- E->sphere.dircos[1][3] = -sin(ro);
- E->sphere.dircos[2][1] = -sin(fo);
- E->sphere.dircos[2][2] = cos(fo);
- E->sphere.dircos[2][3] = 0.0;
- E->sphere.dircos[3][1] = sin(ro)*cos(fo);
- E->sphere.dircos[3][2] = sin(ro)*sin(fo);
- E->sphere.dircos[3][3] = cos(ro);
-
for (j=1;j<=E->sphere.caps_per_proc;j++) {
ii = E->sphere.capid[j];
full_coord_of_cap(E,j,ii);
}
+
+ if (E->control.verbose) {
+ for (lev=E->mesh.levmin;lev<=E->mesh.levmax;lev++) {
+ fprintf(E->fp_out,"output_coordinates before rotation %d \n",lev);
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ for (i=1;i<=E->lmesh.NNO[lev];i++)
+ if(i%E->lmesh.NOZ[lev]==1)
+ fprintf(E->fp_out,"%d %d %g %g %g\n",j,i,E->SX[lev][j][1][i],E->SX[lev][j][2][i],E->SX[lev][j][3][i]);
+ }
+ fflush(E->fp_out);
+ }
+
/* rotate the mesh to avoid two poles on mesh points */
+
+ ro = -0.5*(M_PI/4.0)/E->mesh.elx;
+ fo = 0.0;
+
+ dircos[1][1] = cos(ro)*cos(fo);
+ dircos[1][2] = cos(ro)*sin(fo);
+ dircos[1][3] = -sin(ro);
+ dircos[2][1] = -sin(fo);
+ dircos[2][2] = cos(fo);
+ dircos[2][3] = 0.0;
+ dircos[3][1] = sin(ro)*cos(fo);
+ dircos[3][2] = sin(ro)*sin(fo);
+ dircos[3][3] = cos(ro);
+
for (j=1;j<=E->sphere.caps_per_proc;j++) {
ii = E->sphere.capid[j];
- rotate_mesh(E,j,ii);
+ full_rotate_mesh(E,dircos,j,ii);
}
+ if (E->control.verbose) {
+ for (lev=E->mesh.levmin;lev<=E->mesh.levmax;lev++) {
+ fprintf(E->fp_out,"output_coordinates after rotation %d \n",lev);
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ for (i=1;i<=E->lmesh.NNO[lev];i++)
+ if(i%E->lmesh.NOZ[lev]==1)
+ fprintf(E->fp_out,"%d %d %g %g %g\n",j,i,E->SX[lev][j][1][i],E->SX[lev][j][2][i],E->SX[lev][j][3][i]);
+ }
+ fflush(E->fp_out);
+ }
+
compute_angle_surf_area (E); /* used for interpolation */
for (lev=E->mesh.levmin;lev<=E->mesh.levmax;lev++)
@@ -243,22 +301,6 @@
E->SinCos[lev][j][2][i] = cos(E->SX[lev][j][1][i]);
E->SinCos[lev][j][3][i] = cos(E->SX[lev][j][2][i]);
}
-
- /*
-if (E->control.verbose) {
- for (lev=E->mesh.levmin;lev<=E->mesh.levmax;lev++) {
- fprintf(E->fp_out,"output_coordinates after rotation %d \n",lev);
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- for (i=1;i<=E->lmesh.NNO[lev];i++)
- if(i%E->lmesh.NOZ[lev]==1)
- fprintf(E->fp_out,"%d %d %g %g %g\n",j,i,E->SX[lev][j][1][i],E->SX[lev][j][2][i],E->SX[lev][j][3][i]);
- }
- fflush(E->fp_out);
-}
- */
-
-
-
return;
}
Modified: mc/3D/CitcomS/trunk/lib/Sphere_harmonics.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Sphere_harmonics.c 2007-12-20 20:27:41 UTC (rev 8962)
+++ mc/3D/CitcomS/trunk/lib/Sphere_harmonics.c 2007-12-20 20:30:13 UTC (rev 8963)
@@ -141,195 +141,3 @@
return;
}
-/* ================================================
- compute angle and area
- ================================================*/
-
-void compute_angle_surf_area (E)
- struct All_variables *E;
-{
-
- int es,el,m,i,j,ii,ia[5],lev;
- double aa,y1[4],y2[4],angle[6],xx[4][5],area_sphere_cap();
- void get_angle_sphere_cap();
- void parallel_process_termination();
-
- for (m=1;m<=E->sphere.caps_per_proc;m++) {
- ia[1] = 1;
- ia[2] = E->lmesh.noz*E->lmesh.nox-E->lmesh.noz+1;
- ia[3] = E->lmesh.nno-E->lmesh.noz+1;
- ia[4] = ia[3]-E->lmesh.noz*(E->lmesh.nox-1);
-
- for (i=1;i<=4;i++) {
- xx[1][i] = E->x[m][1][ia[i]]/E->sx[m][3][ia[1]];
- xx[2][i] = E->x[m][2][ia[i]]/E->sx[m][3][ia[1]];
- xx[3][i] = E->x[m][3][ia[i]]/E->sx[m][3][ia[1]];
- }
-
- get_angle_sphere_cap(xx,angle);
-
- for (i=1;i<=4;i++) /* angle1: bet 1 & 2; angle2: bet 2 & 3 ..*/
- E->sphere.angle[m][i] = angle[i];
-
- E->sphere.area[m] = area_sphere_cap(angle);
-
- for (lev=E->mesh.levmax;lev>=E->mesh.levmin;lev--)
- for (es=1;es<=E->lmesh.SNEL[lev];es++) {
- el = (es-1)*E->lmesh.ELZ[lev]+1;
- for (i=1;i<=4;i++)
- ia[i] = E->IEN[lev][m][el].node[i];
-
- for (i=1;i<=4;i++) {
- xx[1][i] = E->X[lev][m][1][ia[i]]/E->SX[lev][m][3][ia[1]];
- xx[2][i] = E->X[lev][m][2][ia[i]]/E->SX[lev][m][3][ia[1]];
- xx[3][i] = E->X[lev][m][3][ia[i]]/E->SX[lev][m][3][ia[1]];
- }
-
- get_angle_sphere_cap(xx,angle);
-
- for (i=1;i<=4;i++) /* angle1: bet 1 & 2; angle2: bet 2 & 3 ..*/
- E->sphere.angle1[lev][m][i][es] = angle[i];
-
- E->sphere.area1[lev][m][es] = area_sphere_cap(angle);
-
-/* fprintf(E->fp_out,"lev%d %d %.6e %.6e %.6e %.6e %.6e\n",lev,es,angle[1],angle[2],angle[3],angle[4],E->sphere.area1[lev][m][es]); */
-
- } /* end for lev and es */
-
- } /* end for m */
-
- return;
-}
-
-/* ================================================
- area of spherical rectangle
- ================================================ */
-double area_sphere_cap(angle)
- double angle[6];
-{
-
- double area,a,b,c;
- double area_of_sphere_triag();
-
- a = angle[1];
- b = angle[2];
- c = angle[5];
- area = area_of_sphere_triag(a,b,c);
-
- a = angle[3];
- b = angle[4];
- c = angle[5];
- area += area_of_sphere_triag(a,b,c);
-
- return (area);
-}
-
-/* ================================================
- area of spherical triangle
- ================================================ */
-double area_of_sphere_triag(a,b,c)
- double a,b,c;
-{
-
- double ss,ak,aa,bb,cc,area;
- const double e_16 = 1.0e-16;
- const double two = 2.0;
- const double pt5 = 0.5;
-
- ss = (a+b+c)*pt5;
- area=0.0;
- a = sin(ss-a);
- b = sin(ss-b);
- c = sin(ss-c);
- ak = a*b*c/sin(ss); /* sin(ss-a)*sin(ss-b)*sin(ss-c)/sin(ss) */
- if(ak<e_16) return (area);
- ak = sqrt(ak);
- aa = two*atan(ak/a);
- bb = two*atan(ak/b);
- cc = two*atan(ak/c);
- area = aa+bb+cc-M_PI;
-
- return (area);
-}
-
-/* =====================================================================
- get the area for given five points (4 nodes for a rectangle and one test node)
- angle [i]: angle bet test node and node i of the rectangle
- angle1[i]: angle bet nodes i and i+1 of the rectangle
- ====================================================================== */
-double area_of_5points(E,lev,m,el,x,ne)
- struct All_variables *E;
- int lev,m,el,ne;
- double x[4];
-{
- int i,es,ia[5];
- double area1,get_angle(),area_of_sphere_triag();
- double xx[4],angle[5],angle1[5];
-
- for (i=1;i<=4;i++)
- ia[i] = E->IEN[lev][m][el].node[i];
-
- es = (el-1)/E->lmesh.ELZ[lev]+1;
-
- for (i=1;i<=4;i++) {
- xx[1] = E->X[lev][m][1][ia[i]]/E->SX[lev][m][3][ia[1]];
- xx[2] = E->X[lev][m][2][ia[i]]/E->SX[lev][m][3][ia[1]];
- xx[3] = E->X[lev][m][3][ia[i]]/E->SX[lev][m][3][ia[1]];
- angle[i] = get_angle(x,xx); /* get angle bet (i,j) and other four*/
- angle1[i]= E->sphere.angle1[lev][m][i][es];
- }
-
- area1 = area_of_sphere_triag(angle[1],angle[2],angle1[1])
- + area_of_sphere_triag(angle[2],angle[3],angle1[2])
- + area_of_sphere_triag(angle[3],angle[4],angle1[3])
- + area_of_sphere_triag(angle[4],angle[1],angle1[4]);
-
- return (area1);
-}
-
-/* ================================
- get the angle for given four points spherical rectangle
- ================================= */
-
-void get_angle_sphere_cap(xx,angle)
- double xx[4][5],angle[6];
-{
-
- int i,j,ii;
- double y1[4],y2[4],get_angle();;
-
- for (i=1;i<=4;i++) { /* angle1: bet 1 & 2; angle2: bet 2 & 3 ..*/
- for (j=1;j<=3;j++) {
- ii=(i==4)?1:(i+1);
- y1[j] = xx[j][i];
- y2[j] = xx[j][ii];
- }
- angle[i] = get_angle(y1,y2);
- }
-
- for (j=1;j<=3;j++) {
- y1[j] = xx[j][1];
- y2[j] = xx[j][3];
- }
-
- angle[5] = get_angle(y1,y2); /* angle5 for betw 1 and 3: diagonal */
- return;
-}
-
-/* ================================
- get the angle for given two points
- ================================= */
-double get_angle(x,xx)
- double x[4],xx[4];
-{
- double dist,angle;
- const double pt5 = 0.5;
- const double two = 2.0;
-
- dist=sqrt( (x[1]-xx[1])*(x[1]-xx[1])
- + (x[2]-xx[2])*(x[2]-xx[2])
- + (x[3]-xx[3])*(x[3]-xx[3]) )*pt5;
- angle = asin(dist)*two;
-
- return (angle);
-}
Modified: mc/3D/CitcomS/trunk/lib/Sphere_util.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Sphere_util.c 2007-12-20 20:27:41 UTC (rev 8962)
+++ mc/3D/CitcomS/trunk/lib/Sphere_util.c 2007-12-20 20:30:13 UTC (rev 8963)
@@ -60,37 +60,195 @@
return;
}
+/* ================================================
+ compute angle and area
+ ================================================*/
-/* =================================================
- rotate the mesh
- =================================================*/
-void rotate_mesh(E,m,icap)
- struct All_variables *E;
- int icap,m;
- {
+void compute_angle_surf_area (E)
+ struct All_variables *E;
+{
- int i,lev;
- double t[4],myatan();
+ int es,el,m,i,j,ii,ia[5],lev;
+ double aa,y1[4],y2[4],angle[6],xx[4][5],area_sphere_cap();
+ void get_angle_sphere_cap();
+ void parallel_process_termination();
- for (lev=E->mesh.levmin;lev<=E->mesh.levmax;lev++) {
- for (i=1;i<=E->lmesh.NNO[lev];i++) {
- t[0] = E->X[lev][m][1][i]*E->sphere.dircos[1][1]+
- E->X[lev][m][2][i]*E->sphere.dircos[1][2]+
- E->X[lev][m][3][i]*E->sphere.dircos[1][3];
- t[1] = E->X[lev][m][1][i]*E->sphere.dircos[2][1]+
- E->X[lev][m][2][i]*E->sphere.dircos[2][2]+
- E->X[lev][m][3][i]*E->sphere.dircos[2][3];
- t[2] = E->X[lev][m][1][i]*E->sphere.dircos[3][1]+
- E->X[lev][m][2][i]*E->sphere.dircos[3][2]+
- E->X[lev][m][3][i]*E->sphere.dircos[3][3];
+ for (m=1;m<=E->sphere.caps_per_proc;m++) {
+ ia[1] = 1;
+ ia[2] = E->lmesh.noz*E->lmesh.nox-E->lmesh.noz+1;
+ ia[3] = E->lmesh.nno-E->lmesh.noz+1;
+ ia[4] = ia[3]-E->lmesh.noz*(E->lmesh.nox-1);
- E->X[lev][m][1][i] = t[0];
- E->X[lev][m][2][i] = t[1];
- E->X[lev][m][3][i] = t[2];
- E->SX[lev][m][1][i] = acos(t[2]/E->SX[lev][m][3][i]);
- E->SX[lev][m][2][i] = myatan(t[1],t[0]);
- }
- } /* lev */
+ for (i=1;i<=4;i++) {
+ xx[1][i] = E->x[m][1][ia[i]]/E->sx[m][3][ia[1]];
+ xx[2][i] = E->x[m][2][ia[i]]/E->sx[m][3][ia[1]];
+ xx[3][i] = E->x[m][3][ia[i]]/E->sx[m][3][ia[1]];
+ }
- return;
- }
+ get_angle_sphere_cap(xx,angle);
+
+ for (i=1;i<=4;i++) /* angle1: bet 1 & 2; angle2: bet 2 & 3 ..*/
+ E->sphere.angle[m][i] = angle[i];
+
+ E->sphere.area[m] = area_sphere_cap(angle);
+
+ for (lev=E->mesh.levmax;lev>=E->mesh.levmin;lev--)
+ for (es=1;es<=E->lmesh.SNEL[lev];es++) {
+ el = (es-1)*E->lmesh.ELZ[lev]+1;
+ for (i=1;i<=4;i++)
+ ia[i] = E->IEN[lev][m][el].node[i];
+
+ for (i=1;i<=4;i++) {
+ xx[1][i] = E->X[lev][m][1][ia[i]]/E->SX[lev][m][3][ia[1]];
+ xx[2][i] = E->X[lev][m][2][ia[i]]/E->SX[lev][m][3][ia[1]];
+ xx[3][i] = E->X[lev][m][3][ia[i]]/E->SX[lev][m][3][ia[1]];
+ }
+
+ get_angle_sphere_cap(xx,angle);
+
+ for (i=1;i<=4;i++) /* angle1: bet 1 & 2; angle2: bet 2 & 3 ..*/
+ E->sphere.angle1[lev][m][i][es] = angle[i];
+
+ E->sphere.area1[lev][m][es] = area_sphere_cap(angle);
+
+/* fprintf(E->fp_out,"lev%d %d %.6e %.6e %.6e %.6e %.6e\n",lev,es,angle[1],angle[2],angle[3],angle[4],E->sphere.area1[lev][m][es]); */
+
+ } /* end for lev and es */
+
+ } /* end for m */
+
+ return;
+}
+
+/* ================================================
+ area of spherical rectangle
+ ================================================ */
+double area_sphere_cap(angle)
+ double angle[6];
+{
+
+ double area,a,b,c;
+ double area_of_sphere_triag();
+
+ a = angle[1];
+ b = angle[2];
+ c = angle[5];
+ area = area_of_sphere_triag(a,b,c);
+
+ a = angle[3];
+ b = angle[4];
+ c = angle[5];
+ area += area_of_sphere_triag(a,b,c);
+
+ return (area);
+}
+
+/* ================================================
+ area of spherical triangle
+ ================================================ */
+double area_of_sphere_triag(a,b,c)
+ double a,b,c;
+{
+
+ double ss,ak,aa,bb,cc,area;
+ const double e_16 = 1.0e-16;
+ const double two = 2.0;
+ const double pt5 = 0.5;
+
+ ss = (a+b+c)*pt5;
+ area=0.0;
+ a = sin(ss-a);
+ b = sin(ss-b);
+ c = sin(ss-c);
+ ak = a*b*c/sin(ss); /* sin(ss-a)*sin(ss-b)*sin(ss-c)/sin(ss) */
+ if(ak<e_16) return (area);
+ ak = sqrt(ak);
+ aa = two*atan(ak/a);
+ bb = two*atan(ak/b);
+ cc = two*atan(ak/c);
+ area = aa+bb+cc-M_PI;
+
+ return (area);
+}
+
+/* =====================================================================
+ get the area for given five points (4 nodes for a rectangle and one test node)
+ angle [i]: angle bet test node and node i of the rectangle
+ angle1[i]: angle bet nodes i and i+1 of the rectangle
+ ====================================================================== */
+double area_of_5points(E,lev,m,el,x,ne)
+ struct All_variables *E;
+ int lev,m,el,ne;
+ double x[4];
+{
+ int i,es,ia[5];
+ double area1,get_angle(),area_of_sphere_triag();
+ double xx[4],angle[5],angle1[5];
+
+ for (i=1;i<=4;i++)
+ ia[i] = E->IEN[lev][m][el].node[i];
+
+ es = (el-1)/E->lmesh.ELZ[lev]+1;
+
+ for (i=1;i<=4;i++) {
+ xx[1] = E->X[lev][m][1][ia[i]]/E->SX[lev][m][3][ia[1]];
+ xx[2] = E->X[lev][m][2][ia[i]]/E->SX[lev][m][3][ia[1]];
+ xx[3] = E->X[lev][m][3][ia[i]]/E->SX[lev][m][3][ia[1]];
+ angle[i] = get_angle(x,xx); /* get angle bet (i,j) and other four*/
+ angle1[i]= E->sphere.angle1[lev][m][i][es];
+ }
+
+ area1 = area_of_sphere_triag(angle[1],angle[2],angle1[1])
+ + area_of_sphere_triag(angle[2],angle[3],angle1[2])
+ + area_of_sphere_triag(angle[3],angle[4],angle1[3])
+ + area_of_sphere_triag(angle[4],angle[1],angle1[4]);
+
+ return (area1);
+}
+
+/* ================================
+ get the angle for given four points spherical rectangle
+ ================================= */
+
+void get_angle_sphere_cap(xx,angle)
+ double xx[4][5],angle[6];
+{
+
+ int i,j,ii;
+ double y1[4],y2[4],get_angle();;
+
+ for (i=1;i<=4;i++) { /* angle1: bet 1 & 2; angle2: bet 2 & 3 ..*/
+ for (j=1;j<=3;j++) {
+ ii=(i==4)?1:(i+1);
+ y1[j] = xx[j][i];
+ y2[j] = xx[j][ii];
+ }
+ angle[i] = get_angle(y1,y2);
+ }
+
+ for (j=1;j<=3;j++) {
+ y1[j] = xx[j][1];
+ y2[j] = xx[j][3];
+ }
+
+ angle[5] = get_angle(y1,y2); /* angle5 for betw 1 and 3: diagonal */
+ return;
+}
+
+/* ================================
+ get the angle for given two points
+ ================================= */
+double get_angle(x,xx)
+ double x[4],xx[4];
+{
+ double dist,angle;
+ const double pt5 = 0.5;
+ const double two = 2.0;
+
+ dist=sqrt( (x[1]-xx[1])*(x[1]-xx[1])
+ + (x[2]-xx[2])*(x[2]-xx[2])
+ + (x[3]-xx[3])*(x[3]-xx[3]) )*pt5;
+ angle = asin(dist)*two;
+
+ return (angle);
+}
Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2007-12-20 20:27:41 UTC (rev 8962)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2007-12-20 20:30:13 UTC (rev 8963)
@@ -276,7 +276,6 @@
double *area1[MAX_LEVELS][NCS];
double *angle1[MAX_LEVELS][NCS][5];
- double dircos[4][4];
double *R[MAX_LEVELS];
double ro,ri;
struct CAP cap[NCS];
Added: mc/3D/CitcomS/trunk/tests/gnomonic.c
===================================================================
--- mc/3D/CitcomS/trunk/tests/gnomonic.c 2007-12-20 20:27:41 UTC (rev 8962)
+++ mc/3D/CitcomS/trunk/tests/gnomonic.c 2007-12-20 20:30:13 UTC (rev 8963)
@@ -0,0 +1,71 @@
+/*
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ *
+ *<LicenseText>
+ *
+ * CitcomS by Louis Moresi, Shijie Zhong, Lijie Han, Eh Tan,
+ * Clint Conrad, Michael Gurnis, and Eun-seo Choi.
+ * Copyright (C) 1994-2005, California Institute of Technology.
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ *
+ *</LicenseText>
+ *
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ */
+
+
+/* Compile by: mpicc gnomonic.c -lCitcomS -L../lib -lz */
+
+#include <math.h>
+#include <stdio.h>
+
+void spherical_to_uv2(double center[2], int len,
+ double *theta, double *phi,
+ double *u, double *v);
+void uv_to_spherical(double center[2], int len,
+ double *u, double *v,
+ double *theta, double *phi);
+
+/* test for gnomonic projection and its inverse */
+
+
+int main(int argc, char **argv)
+{
+ #define len 6
+ int i;
+
+ double u[len], v[len];
+ double center[2] = {M_PI / 2, 0};
+ double theta[len] = {0.1, 0.2, 0.3, 0.3, center[0], M_PI/4};
+ double phi[len] = {0.1, 0.1, 0.3, 0.4, center[1], M_PI/4};
+
+ spherical_to_uv2(center, len, theta, phi, u, v);
+
+ for(i=0; i<len; i++) {
+ fprintf(stderr, "(%.15e, %.15e) -> (%.15e, %.15e)\n",
+ theta[i], phi[i], u[i], v[i]);
+ }
+ fprintf(stderr, "\n\n");
+
+ uv_to_spherical(center, len, u, v, theta, phi);
+
+ for(i=0; i<len; i++) {
+ fprintf(stderr, "(%.15e, %.15e) <- (%.15e, %.15e)\n",
+ theta[i], phi[i], u[i], v[i]);
+ }
+
+ return 0;
+}
More information about the cig-commits
mailing list