/*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*
*<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>
*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/
/* Functions relating to the building and use of mesh locations ... */
#include <math.h>
#include <sys/types.h>
#include "element_definitions.h"
#include "global_defs.h"
/**************************************************************/
/* 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 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 icap )
{
int i, j, k, lev, temp, elx, ely;
int node, snode, ns, step;
int lelx, lely, lnox, lnoy;
int lvnox, lvnoy, lvnoz;
int ok;
double x[5], y[5], z[5], center[3], referencep[2];
double xx[5], yy[5];
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, efac2,rfac;
double a, b;
double offset;
double myatan();
void even_divide_arc12();
temp = max(E->mesh.noy, E->mesh.nox);
theta0 = (double *)malloc((temp+1)*sizeof(double));
fi0 = (double *)malloc((temp+1)*sizeof(double));
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));
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.nox;
px = malloc((temp+1)*sizeof(double));
py = malloc((temp+1)*sizeof(double));
qx = malloc((temp+1)*sizeof(double));
qy = malloc((temp+1)*sizeof(double));
/* define the cap corners */
/* adjust the corner coordinates so that the size (surface area) of
each cap is about the same. */
offset = 9.736/180.0*M_PI;
for (i=1;i<=4;i++) {
E->sphere.cap[(i-1)*3+1].theta[1] = 0.0;
E->sphere.cap[(i-1)*3+1].theta[2] = M_PI/4.0+offset;
E->sphere.cap[(i-1)*3+1].theta[3] = M_PI/2.0;
E->sphere.cap[(i-1)*3+1].theta[4] = M_PI/4.0+offset;
E->sphere.cap[(i-1)*3+1].fi[1] = 0.0;
E->sphere.cap[(i-1)*3+1].fi[2] = (i-1)*M_PI/2.0;
E->sphere.cap[(i-1)*3+1].fi[3] = (i-1)*M_PI/2.0 + M_PI/4.0;
E->sphere.cap[(i-1)*3+1].fi[4] = i*M_PI/2.0;
E->sphere.cap[(i-1)*3+2].theta[1] = M_PI/4.0+offset;
E->sphere.cap[(i-1)*3+2].theta[2] = M_PI/2.0;
E->sphere.cap[(i-1)*3+2].theta[3] = 3*M_PI/4.0-offset;
E->sphere.cap[(i-1)*3+2].theta[4] = M_PI/2.0;
E->sphere.cap[(i-1)*3+2].fi[1] = i*M_PI/2.0;
E->sphere.cap[(i-1)*3+2].fi[2] = i*M_PI/2.0 - M_PI/4.0;
E->sphere.cap[(i-1)*3+2].fi[3] = i*M_PI/2.0;
E->sphere.cap[(i-1)*3+2].fi[4] = i*M_PI/2.0 + M_PI/4.0;
}
for (i=1;i<=4;i++) {
j = (i-1)*3;
if (i==1) j=12;
E->sphere.cap[j].theta[1] = M_PI/2.0;
E->sphere.cap[j].theta[2] = 3*M_PI/4.0-offset;
E->sphere.cap[j].theta[3] = M_PI;
E->sphere.cap[j].theta[4] = 3*M_PI/4.0-offset;
E->sphere.cap[j].fi[1] = (i-1)*M_PI/2.0 + M_PI/4.0;
E->sphere.cap[j].fi[2] = (i-1)*M_PI/2.0;
E->sphere.cap[j].fi[3] = 0.0;
E->sphere.cap[j].fi[4] = i*M_PI/2.0;
}
/* 4 corners of the cap in Cartesian coordinates */
/* the order of corners is: */
/* 1 - 4 */
/* | | */
/* 2 - 3 */
center[0] = center[1] = center[2] = 0;
#ifdef ALLOW_ELLIPTICAL
for (i=1;i<=4;i++) { /* works for both elliptical and spherical */
x[i] = E->data.ra * sin(E->sphere.cap[icap].theta[i])*cos(E->sphere.cap[icap].fi[i]);
y[i] = E->data.ra * sin(E->sphere.cap[icap].theta[i])*sin(E->sphere.cap[icap].fi[i]);
z[i] = E->data.rc * cos(E->sphere.cap[icap].theta[i]);
center[0] += x[i];
center[1] += y[i];
center[2] += z[i];
}
#else
/* only spherical */
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]);
center[0] += x[i];
center[1] += y[i];
center[2] += z[i];
}
#endif
center[0] *= 0.25;
center[1] *= 0.25;
center[2] *= 0.25;
/* use the center as the reference point for gnomonic projection */
referencep[0] = acos(center[2] /
sqrt(center[0]*center[0] +
center[1]*center[1] +
center[2]*center[2]));;
referencep[1] = myatan(center[1], center[0]);
lev = E->mesh.levmax;
/* # of elements/nodes per cap */
elx = E->lmesh.ELX[lev]*E->parallel.nprocx;
ely = E->lmesh.ELY[lev]*E->parallel.nprocy;
/* # of elements/nodes per processor */
lelx = E->lmesh.ELX[lev];
lely = E->lmesh.ELY[lev];
lnox = lelx+1;
lnoy = lely+1;
/* evenly divide arc linking corner 1 and 2 */
even_divide_arc12(elx,x[1],y[1],z[1],x[2],y[2],z[2],theta0,fi0);
/* 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];
}
/* evenly divide arc linking corner 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];
}
/* evenly divide arc linking corner 1 and 4 */
even_divide_arc12(ely,x[1],y[1],z[1],x[4],y[4],z[4],theta0,fi0);
/* 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 corner 2 and 3 */
even_divide_arc12(ely,x[2],y[2],z[2],x[3],y[3],z[3],theta0,fi0);
/* 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];
}
/* 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 */
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);
snode = 0;
for(k=0; k<lnoy; k++) {
xx[2] = u3[k];
yy[2] = v3[k];
xx[3] = u4[k];
yy[3] = v4[k];
for(j=0; j<lnox; j++) {
xx[0] = u1[j];
yy[0] = v1[j];
xx[1] = u2[j];
yy[1] = v2[j];
ok = find_intersection(xx, yy, &a, &b);
if(!ok) {
fprintf(stderr, "Error(Full_coord_of_cap): cannot find intersection point! rank=%d, nx=%d, ny=%d\n", E->parallel.me, j, k);
fprintf(stderr, "L1: (%g, %g)-(%g, %g) L2 (%g, %g)-(%g, %g)\n",
xx[0],yy[0],xx[1],yy[1],xx[2],yy[2],xx[3],yy[3]);
exit(10);
}
px[snode] = a;
py[snode] = b;
snode++;
}
}
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];
}
}
#ifdef ALLOW_ELLIPTICAL
/* both spherical and elliptical */
efac2 = E->data.ellipticity*(2.0 - E->data.ellipticity)/
((1.- E->data.ellipticity)*(1.-E->data.ellipticity));
for (lev=E->mesh.levmax, step=1; lev>=E->mesh.levmin; lev--, step*=2) {
/* store the node location in spherical and cartesian coordinates */
lvnox = E->lmesh.NOX[lev];
lvnoy = E->lmesh.NOY[lev];
lvnoz = E->lmesh.NOZ[lev];
node = 1;
for (k=0; k<lvnoy; k++) {
for (j=0, ns=step*lnoy*k; j<lvnox; j++, ns+=step) {
theta = qx[ns];
fi = qy[ns];
cost = cos(theta);
rfac = E->data.ra*1./sqrt(1.0+efac2*cost*cost);
sint = sin(theta);
cosf = cos(fi);
sinf = sin(fi);
for (i=1; i<=lvnoz; i++) {
/* theta,fi,and r coordinates */
E->SX[lev][1][node] = theta;
E->SX[lev][2][node] = fi;
E->SX[lev][3][node] = rfac * E->sphere.R[lev][i];
/* x,y,and z oordinates */
E->X[lev][1][node] = E->data.ra * E->sphere.R[lev][i]*sint*cosf;
E->X[lev][2][node] = E->data.ra * E->sphere.R[lev][i]*sint*sinf;
E->X[lev][3][node] = E->data.rc * E->sphere.R[lev][i]*cost;
node++;
}
}
}
} /* end for lev */
#else
/* spherical */
for (lev=E->mesh.levmax, step=1; lev>=E->mesh.levmin; lev--, step*=2) {
/* store the node location in spherical and cartesian coordinates */
lvnox = E->lmesh.NOX[lev];
lvnoy = E->lmesh.NOY[lev];
lvnoz = E->lmesh.NOZ[lev];
node = 1;
for (k=0; k<lvnoy; k++) {
for (j=0, ns=step*lnoy*k; j<lvnox; j++, ns+=step) {
theta = qx[ns];
fi = qy[ns];
cost = cos(theta);
sint = sin(theta);
cosf = cos(fi);
sinf = sin(fi);
for (i=1; i<=lvnoz; i++) {
/* theta,fi,and r coordinates */
E->SX[lev][1][node] = theta;
E->SX[lev][2][node] = fi;
E->SX[lev][3][node] = E->sphere.R[lev][i];
/* x,y,and z oordinates */
E->X[lev][1][node] = E->sphere.R[lev][i]*sint*cosf;
E->X[lev][2][node] = E->sphere.R[lev][i]*sint*sinf;
E->X[lev][3][node] = E->sphere.R[lev][i]*cost;
node++;
}
}
}
} /* end for lev */
#endif
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 );
}