https://github.com/geodynamics/citcoms
Revision ebff88cbd36ec28feb0294775d6e10e70fccaf21 authored by Eh Tan on 16 March 2009, 22:52:08 UTC, committed by Eh Tan on 16 March 2009, 22:52:08 UTC
1 parent e1f9320
Tip revision: ebff88cbd36ec28feb0294775d6e10e70fccaf21 authored by Eh Tan on 16 March 2009, 22:52:08 UTC
add test cases for the stokes solver under bousinessq approximation and truncated anelastic liquid approximation
add test cases for the stokes solver under bousinessq approximation and truncated anelastic liquid approximation
Tip revision: ebff88c
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 ...