https://github.com/geodynamics/citcoms
Revision df821c06952128812763ec5b7a944cd178ca5ed3 authored by Thorsten Becker on 08 December 2010, 18:43:58 UTC, committed by Thorsten Becker on 08 December 2010, 18:43:58 UTC
GMT/ggrd compile before).
1 parent cf7495e
Tip revision: df821c06952128812763ec5b7a944cd178ca5ed3 authored by Thorsten Becker on 08 December 2010, 18:43:58 UTC
Added definition of TRUE (1) and FALSE (0) in case undefined (was defined for
Added definition of TRUE (1) and FALSE (0) in case undefined (was defined for
Tip revision: df821c0
Sphere_util.c
/*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*
*<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 <math.h>
#include "global_defs.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(elx,x1,y1,z1,x2,y2,z2,theta,fi)
double x1,y1,z1,x2,y2,z2,*theta,*fi;
int elx;
{
double dx,dy,dz,xx,yy,zz,myatan();
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 (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);
}
Computing file changes ...