https://github.com/geodynamics/citcoms
Revision 36f5b324acac508e900f74fcb6f9c7ae35fc8cbc authored by Leif Strand on 29 April 2009, 20:47:56 UTC, committed by Leif Strand on 29 April 2009, 20:47:56 UTC
1 parent 5ceef5d
Tip revision: 36f5b324acac508e900f74fcb6f9c7ae35fc8cbc authored by Leif Strand on 29 April 2009, 20:47:56 UTC
Merged r14275:14351 from trunk.
Merged r14275:14351 from trunk.
Tip revision: 36f5b32
Sphere_util.cc
/*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*
*<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>
*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/
/* Common functions relating to the building and use of mesh locations ... */
#include "sphere_util.h"
#include <math.h>
#include "global_defs.h"
#include "pan_problem_misc_functions.h"
/* =================================================
this routine evenly divides the arc between points
1 and 2 in a great cicle. The word "evenly" means
anglewise evenly.
=================================================*/
void even_divide_arc12(
int elx,
double x1, double y1, double z1, double x2, double y2, double z2, double *theta, double *fi
)
{
double dx,dy,dz,xx,yy,zz;
int j, nox;
nox = elx+1;
dx = (x2 - x1)/elx;
dy = (y2 - y1)/elx;
dz = (z2 - z1)/elx;
for (j=1;j<=nox;j++) {
xx = x1 + dx*(j-1) + 5.0e-32;
yy = y1 + dy*(j-1);
zz = z1 + dz*(j-1);
theta[j] = acos(zz/sqrt(xx*xx+yy*yy+zz*zz));
fi[j] = myatan(yy,xx);
}
return;
}
/* ================================================
compute angle and area
================================================*/
void compute_angle_surf_area(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];
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(double angle[6])
{
double area,a,b,c;
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(double a, double b, double 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(
struct All_variables *E,
int lev, int m, int el,
double x[4],
int ne
)
{
int i,es,ia[5];
double area1;
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(double xx[4][5], double angle[6])
{
int i,j,ii;
double y1[4],y2[4];
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(double x[4], double 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);
}
Computing file changes ...