/* *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ * * * * 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 * * * *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */ /* Common functions relating to the building and use of mesh locations ... */ #include #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(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]; 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(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(akIEN[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); }