[cig-commits] r7813 - mc/3D/CitcomS/trunk/lib
tan2 at geodynamics.org
tan2 at geodynamics.org
Mon Aug 13 15:55:42 PDT 2007
Author: tan2
Date: 2007-08-13 15:55:41 -0700 (Mon, 13 Aug 2007)
New Revision: 7813
Modified:
mc/3D/CitcomS/trunk/lib/Global_operations.c
mc/3D/CitcomS/trunk/lib/Nodal_mesh.c
mc/3D/CitcomS/trunk/lib/Obsolete.c
mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c
mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c
Log:
Minor refactoring, clean-up and moving functions between files
Modified: mc/3D/CitcomS/trunk/lib/Global_operations.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Global_operations.c 2007-08-13 22:04:33 UTC (rev 7812)
+++ mc/3D/CitcomS/trunk/lib/Global_operations.c 2007-08-13 22:55:41 UTC (rev 7813)
@@ -64,9 +64,19 @@
}
return;
- }
+}
+void remove_horiz_ave2(struct All_variables *E, double **X)
+{
+ double *H;
+
+ H = (double *)malloc( (E->lmesh.noz+1)*sizeof(double));
+ remove_horiz_ave(E, X, H, 0);
+ free ((void *) H);
+}
+
+
void return_horiz_ave(E,X,H)
struct All_variables *E;
double **X, *H;
Modified: mc/3D/CitcomS/trunk/lib/Nodal_mesh.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Nodal_mesh.c 2007-08-13 22:04:33 UTC (rev 7812)
+++ mc/3D/CitcomS/trunk/lib/Nodal_mesh.c 2007-08-13 22:55:41 UTC (rev 7813)
@@ -1,6 +1,6 @@
/*
- *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- *
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ *
*<LicenseText>
*
* CitcomS by Louis Moresi, Shijie Zhong, Lijie Han, Eh Tan,
@@ -22,8 +22,8 @@
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*</LicenseText>
- *
- *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ *
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/
/* Functions relating to the building and use of mesh locations ... */
@@ -33,20 +33,157 @@
#include "element_definitions.h"
#include "global_defs.h"
-extern int Emergency_stop;
-/*
-void flogical_mesh_to_real(E,data,level)
+void v_from_vector(E)
struct All_variables *E;
- float *data;
- int level;
+{
+ int m,node;
+ const int nno = E->lmesh.nno;
-{ int i,j,n1,n2;
+ for (m=1;m<=E->sphere.caps_per_proc;m++) {
+ for(node=1;node<=nno;node++) {
+ E->sphere.cap[m].V[1][node] = E->U[m][E->id[m][node].doff[1]];
+ E->sphere.cap[m].V[2][node] = E->U[m][E->id[m][node].doff[2]];
+ E->sphere.cap[m].V[3][node] = E->U[m][E->id[m][node].doff[3]];
+ if (E->node[m][node] & VBX)
+ E->sphere.cap[m].V[1][node] = E->sphere.cap[m].VB[1][node];
+ if (E->node[m][node] & VBY)
+ E->sphere.cap[m].V[2][node] = E->sphere.cap[m].VB[2][node];
+ if (E->node[m][node] & VBZ)
+ E->sphere.cap[m].V[3][node] = E->sphere.cap[m].VB[3][node];
+ }
+ }
- return;
+ return;
}
-*/
+void v_from_vector_pseudo_surf(E)
+ struct All_variables *E;
+{
+ int m,node;
+
+ const int nno = E->lmesh.nno;
+ double sum_V = 0.0, sum_dV = 0.0, rel_error = 0.0, global_max_error = 0.0;
+ double tol_error = 1.0e-03;
+
+ for (m=1;m<=E->sphere.caps_per_proc;m++) {
+ for(node=1;node<=nno;node++) {
+ E->sphere.cap[m].Vprev[1][node] = E->sphere.cap[m].V[1][node];
+ E->sphere.cap[m].Vprev[2][node] = E->sphere.cap[m].V[2][node];
+ E->sphere.cap[m].Vprev[3][node] = E->sphere.cap[m].V[3][node];
+
+ E->sphere.cap[m].V[1][node] = E->U[m][E->id[m][node].doff[1]];
+ E->sphere.cap[m].V[2][node] = E->U[m][E->id[m][node].doff[2]];
+ E->sphere.cap[m].V[3][node] = E->U[m][E->id[m][node].doff[3]];
+ if (E->node[m][node] & VBX)
+ E->sphere.cap[m].V[1][node] = E->sphere.cap[m].VB[1][node];
+ if (E->node[m][node] & VBY)
+ E->sphere.cap[m].V[2][node] = E->sphere.cap[m].VB[2][node];
+ if (E->node[m][node] & VBZ)
+ E->sphere.cap[m].V[3][node] = E->sphere.cap[m].VB[3][node];
+
+ sum_dV += (E->sphere.cap[m].V[1][node] - E->sphere.cap[m].Vprev[1][node])*(E->sphere.cap[m].V[1][node] - E->sphere.cap[m].Vprev[1][node])
+ + (E->sphere.cap[m].V[2][node] - E->sphere.cap[m].Vprev[2][node])*(E->sphere.cap[m].V[2][node] - E->sphere.cap[m].Vprev[2][node])
+ + (E->sphere.cap[m].V[3][node] - E->sphere.cap[m].Vprev[3][node])*(E->sphere.cap[m].V[3][node] - E->sphere.cap[m].Vprev[3][node]);
+ sum_V += E->sphere.cap[m].V[1][node]*E->sphere.cap[m].V[1][node]
+ + E->sphere.cap[m].V[2][node]*E->sphere.cap[m].V[2][node]
+ + E->sphere.cap[m].V[3][node]*E->sphere.cap[m].V[3][node];
+ }
+ rel_error = sqrt(sum_dV)/sqrt(sum_V);
+ MPI_Allreduce(&rel_error,&global_max_error,1,MPI_DOUBLE,MPI_MAX,E->parallel.world);
+ if(global_max_error <= tol_error) E->monitor.stop_topo_loop = 1;
+ if(E->parallel.me==0)
+ fprintf(stderr,"global_max_error=%e stop_topo_loop=%d\n",global_max_error,E->monitor.stop_topo_loop);
+
+ }
+
+ return;
+}
+
+void velo_from_element(E,VV,m,el,sphere_key)
+ struct All_variables *E;
+ float VV[4][9];
+ int el,m,sphere_key;
+{
+
+ int a, node;
+ double sint, cost, sinf, cosf;
+ const int ends=enodes[E->mesh.nsd];
+ const int lev=E->mesh.levmax;
+
+ if (sphere_key)
+ for(a=1;a<=ends;a++) {
+ node = E->ien[m][el].node[a];
+ VV[1][a] = E->sphere.cap[m].V[1][node];
+ VV[2][a] = E->sphere.cap[m].V[2][node];
+ VV[3][a] = E->sphere.cap[m].V[3][node];
+ }
+ else {
+ for(a=1;a<=ends;a++) {
+ node = E->ien[m][el].node[a];
+
+ sint = E->SinCos[lev][m][0][node];
+ sinf = E->SinCos[lev][m][1][node];
+ cost = E->SinCos[lev][m][2][node];
+ cosf = E->SinCos[lev][m][3][node];
+
+ VV[1][a] = E->sphere.cap[m].V[1][node]*cost*cosf
+ - E->sphere.cap[m].V[2][node]*sinf
+ + E->sphere.cap[m].V[3][node]*sint*cosf;
+ VV[2][a] = E->sphere.cap[m].V[1][node]*cost*sinf
+ + E->sphere.cap[m].V[2][node]*cosf
+ + E->sphere.cap[m].V[3][node]*sint*sinf;
+ VV[3][a] = -E->sphere.cap[m].V[1][node]*sint
+ + E->sphere.cap[m].V[3][node]*cost;
+ }
+ }
+ return;
+}
+
+
+void velo_from_element_d(E,VV,m,el,sphere_key)
+ struct All_variables *E;
+ double VV[4][9];
+ int el,m,sphere_key;
+{
+
+ int a, node;
+ double sint, cost, sinf, cosf;
+ const int dims=E->mesh.nsd;
+ const int ends=enodes[E->mesh.nsd];
+ const int nno=E->lmesh.nno;
+ const int lev=E->mesh.levmax;
+
+ if (sphere_key)
+ for(a=1;a<=ends;a++) {
+ node = E->ien[m][el].node[a];
+ VV[1][a] = E->sphere.cap[m].V[1][node];
+ VV[2][a] = E->sphere.cap[m].V[2][node];
+ VV[3][a] = E->sphere.cap[m].V[3][node];
+ }
+ else {
+ for(a=1;a<=ends;a++) {
+ node = E->ien[m][el].node[a];
+
+ sint = E->SinCos[lev][m][0][node];
+ sinf = E->SinCos[lev][m][1][node];
+ cost = E->SinCos[lev][m][2][node];
+ cosf = E->SinCos[lev][m][3][node];
+
+ VV[1][a] = E->sphere.cap[m].V[1][node]*cost*cosf
+ - E->sphere.cap[m].V[2][node]*sinf
+ + E->sphere.cap[m].V[3][node]*sint*cosf;
+ VV[2][a] = E->sphere.cap[m].V[1][node]*cost*sinf
+ + E->sphere.cap[m].V[2][node]*cosf
+ + E->sphere.cap[m].V[3][node]*sint*sinf;
+ VV[3][a] = -E->sphere.cap[m].V[1][node]*sint
+ + E->sphere.cap[m].V[3][node]*cost;
+ }
+ }
+ return;
+}
+
+
void p_to_nodes(E,P,PN,lev)
struct All_variables *E;
double **P;
@@ -58,97 +195,25 @@
for (m=1;m<=E->sphere.caps_per_proc;m++)
for(node=1;node<=E->lmesh.NNO[lev];node++)
PN[m][node] = 0.0;
-
+
for (m=1;m<=E->sphere.caps_per_proc;m++)
for(element=1;element<=E->lmesh.NEL[lev];element++)
for(j=1;j<=enodes[E->mesh.nsd];j++) {
node = E->IEN[lev][m][element].node[j];
- PN[m][node] += P[m][element] * E->TWW[lev][m][element].node[j] ;
+ PN[m][node] += P[m][element] * E->TWW[lev][m][element].node[j] ;
}
-
+
(E->exchange_node_f)(E,PN,lev);
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(node=1;node<=E->lmesh.NNO[lev];node++)
PN[m][node] *= E->MASS[lev][m][node];
- return;
+ return;
}
-/*
-void p_to_centres(E,PN,P,lev)
- struct All_variables *E;
- float **PN;
- double **P;
- int lev;
-{ int p,element,node,j,m;
- double weight;
-
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(p=1;p<=E->lmesh.NEL[lev];p++)
- P[m][p] = 0.0;
-
- weight=1.0/((double)enodes[E->mesh.nsd]) ;
-
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(p=1;p<=E->lmesh.NEL[lev];p++)
- for(j=1;j<=enodes[E->mesh.nsd];j++)
- P[m][p] += PN[m][E->IEN[lev][m][p].node[j]] * weight;
-
- return;
- }
-*/
-
-/*
-void v_to_intpts(E,VN,VE,lev)
- struct All_variables *E;
- float **VN,**VE;
- int lev;
- {
-
- int m,e,i,j,k;
- const int nsd=E->mesh.nsd;
- const int vpts=vpoints[nsd];
- const int ends=enodes[nsd];
-
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(e=1;e<=E->lmesh.NEL[lev];e++)
- for(i=1;i<=vpts;i++) {
- VE[m][(e-1)*vpts + i] = 0.0;
- for(j=1;j<=ends;j++)
- VE[m][(e-1)*vpts + i] += VN[m][E->IEN[lev][m][e].node[j]]*E->N.vpt[GNVINDEX(j,i)];
- }
-
- return;
- }
-*/
-
-/*
-void visc_to_intpts(E,VN,VE,lev)
- struct All_variables *E;
- float **VN,**VE;
- int lev;
- {
-
- int m,e,i,j,k;
- const int nsd=E->mesh.nsd;
- const int vpts=vpoints[nsd];
- const int ends=enodes[nsd];
-
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(e=1;e<=E->lmesh.NEL[lev];e++)
- for(i=1;i<=vpts;i++) {
- VE[m][(e-1)*vpts + i] = 0.0;
- for(j=1;j<=ends;j++)
- VE[m][(e-1)*vpts + i] += log(VN[m][E->IEN[lev][m][e].node[j]]) * E->N.vpt[GNVINDEX(j,i)];
- VE[m][(e-1)*vpts + i] = exp(VE[m][(e-1)*vpts + i]);
- }
-
- }
-*/
-
void visc_from_gint_to_nodes(E,VE,VN,lev)
struct All_variables *E;
float **VE,**VN;
@@ -176,11 +241,11 @@
VN[m][n] += E->TWW[lev][m][e].node[j] * temp_visc;
}
}
-
+
(E->exchange_node_f)(E,VN,lev);
for(m=1;m<=E->sphere.caps_per_proc;m++)
- for(n=1;n<=E->lmesh.NNO[lev];n++)
+ for(n=1;n<=E->lmesh.NNO[lev];n++)
VN[m][n] *= E->MASS[lev][m][n];
return;
@@ -208,7 +273,7 @@
for(e=1;e<=E->lmesh.NEL[lev];e++)
for(i=1;i<=vpts;i++) {
temp_visc=0.0;
- for(j=1;j<=ends;j++)
+ for(j=1;j<=ends;j++)
temp_visc += E->N.vpt[GNVINDEX(j,i)]*VN[m][E->IEN[lev][m][e].node[j]];
VE[m][(e-1)*vpts+i] = temp_visc;
@@ -241,7 +306,7 @@
VN[m][e] = temp_visc;
}
-
+
return;
}
Modified: mc/3D/CitcomS/trunk/lib/Obsolete.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Obsolete.c 2007-08-13 22:04:33 UTC (rev 7812)
+++ mc/3D/CitcomS/trunk/lib/Obsolete.c 2007-08-13 22:55:41 UTC (rev 7813)
@@ -811,6 +811,91 @@
return;
}
+
+/* ========================================================== */
+/* From Nodal_mesh.c */
+/* =========================================================== */
+
+void flogical_mesh_to_real(E,data,level)
+ struct All_variables *E;
+ float *data;
+ int level;
+
+{ int i,j,n1,n2;
+
+ return;
+}
+
+
+void p_to_centres(E,PN,P,lev)
+ struct All_variables *E;
+ float **PN;
+ double **P;
+ int lev;
+
+{ int p,element,node,j,m;
+ double weight;
+
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ for(p=1;p<=E->lmesh.NEL[lev];p++)
+ P[m][p] = 0.0;
+
+ weight=1.0/((double)enodes[E->mesh.nsd]) ;
+
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ for(p=1;p<=E->lmesh.NEL[lev];p++)
+ for(j=1;j<=enodes[E->mesh.nsd];j++)
+ P[m][p] += PN[m][E->IEN[lev][m][p].node[j]] * weight;
+
+ return;
+ }
+
+
+void v_to_intpts(E,VN,VE,lev)
+ struct All_variables *E;
+ float **VN,**VE;
+ int lev;
+ {
+
+ int m,e,i,j,k;
+ const int nsd=E->mesh.nsd;
+ const int vpts=vpoints[nsd];
+ const int ends=enodes[nsd];
+
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ for(e=1;e<=E->lmesh.NEL[lev];e++)
+ for(i=1;i<=vpts;i++) {
+ VE[m][(e-1)*vpts + i] = 0.0;
+ for(j=1;j<=ends;j++)
+ VE[m][(e-1)*vpts + i] += VN[m][E->IEN[lev][m][e].node[j]]*E->N.vpt[GNVINDEX(j,i)];
+ }
+
+ return;
+ }
+
+
+void visc_to_intpts(E,VN,VE,lev)
+ struct All_variables *E;
+ float **VN,**VE;
+ int lev;
+ {
+
+ int m,e,i,j,k;
+ const int nsd=E->mesh.nsd;
+ const int vpts=vpoints[nsd];
+ const int ends=enodes[nsd];
+
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ for(e=1;e<=E->lmesh.NEL[lev];e++)
+ for(i=1;i<=vpts;i++) {
+ VE[m][(e-1)*vpts + i] = 0.0;
+ for(j=1;j<=ends;j++)
+ VE[m][(e-1)*vpts + i] += log(VN[m][E->IEN[lev][m][e].node[j]]) * E->N.vpt[GNVINDEX(j,i)];
+ VE[m][(e-1)*vpts + i] = exp(VE[m][(e-1)*vpts + i]);
+ }
+
+ }
+
/* version */
/* $Id$ */
Modified: mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c 2007-08-13 22:04:33 UTC (rev 7812)
+++ mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c 2007-08-13 22:55:41 UTC (rev 7813)
@@ -136,11 +136,9 @@
void get_buoyancy(struct All_variables *E, double **buoy)
{
int i,m;
- double temp,temp2,*H;
- void remove_horiz_ave();
+ double temp,temp2;
+ void remove_horiz_ave2(struct All_variables*, double**);
- H = (double *)malloc( (E->lmesh.noz+1)*sizeof(double));
-
temp = E->control.Atemp;
for(m=1;m<=E->sphere.caps_per_proc;m++)
@@ -148,7 +146,7 @@
buoy[m][i] = temp * E->T[m][i];
/* chemical buoyancy */
- if(E->control.tracer &&
+ if(E->control.tracer &&
(E->composition.ichemical_buoyancy)) {
temp2 = E->composition.buoyancy_ratio * temp;
for(m=1;m<=E->sphere.caps_per_proc;m++)
@@ -160,8 +158,7 @@
phase_change_apply_670(E, buoy);
phase_change_apply_cmb(E, buoy);
- remove_horiz_ave(E,buoy,H,0);
- free ((void *) H);
+ remove_horiz_ave2(E,buoy);
return;
}
@@ -394,7 +391,7 @@
return 1.0;
}
-/* convert r,theta,phi system to cartesian, xout[3]
+/* convert r,theta,phi system to cartesian, xout[3]
there's a double version of this in Tracer_setup called
sphere_to_cart
@@ -417,7 +414,7 @@
double ct,cp,st,sp;
-
+
ct=cos(theta);
cp=cos(phi);
st=sin(theta);
@@ -438,7 +435,7 @@
/* given a base from calc_cbase_at_tp, convert a polar vector to
cartesian */
-void convert_pvec_to_cvec(float vr,float vt,
+void convert_pvec_to_cvec(float vr,float vt,
float vp, float *base,
float *cvec)
{
@@ -449,18 +446,18 @@
cvec[i] += base[6+i]* vp;
}
}
-/*
- like malloc, but with test
+/*
+ like malloc, but with test
similar to Malloc1 but I didn't like the int as argument
-
+
*/
void *safe_malloc (size_t size)
{
void *tmp;
-
+
if ((tmp = malloc(size)) == NULL) {
- fprintf(stderr, "safe_malloc: could not allocate memory, %.3f MB\n",
+ fprintf(stderr, "safe_malloc: could not allocate memory, %.3f MB\n",
(float)size/(1024*1024.));
parallel_process_termination();
}
Modified: mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c 2007-08-13 22:04:33 UTC (rev 7812)
+++ mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c 2007-08-13 22:55:41 UTC (rev 7813)
@@ -1,5 +1,5 @@
/*
- *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*
*<LicenseText>
*
@@ -23,8 +23,10 @@
*
*</LicenseText>
*
- *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/
+
+
/* Functions which solve for the velocity and pressure fields using Uzawa-type iteration loop. */
#include <math.h>
@@ -33,481 +35,291 @@
#include "global_defs.h"
#include <stdlib.h>
-extern int Emergency_stop;
-
/* Master loop for pressure and (hence) velocity field */
-
void solve_constrained_flow_iterative(E)
struct All_variables *E;
{
- double *D1;
- double *u;
- double *R,*Bp;
- double residual_ddash;
- double vmag;
- double global_vdot(),global_pdot();
-
float solve_Ahat_p_fhat();
- void assemble_del2_u();
- void assemble_grad_p();
- void assemble_div_u();
void v_from_vector();
void p_to_nodes();
- void strip_bcs_from_residual();
- void velocities_conform_bcs();
- int steps,cycles;
- int i,j,k,doff,vel_cycles_previous,vel_calls_previous;
+ int cycles;
- double time,CPU_time0();
-
- const int npno = E->lmesh.npno;
- const int gnpno = E->mesh.npno;
- const int nno = E->lmesh.nno;
- const int dims = E->mesh.nsd;
- const int neq = E->lmesh.neq;
- const int gneq = E->mesh.neq;
- const int addi_dof = additional_dof[dims];
-
- time=CPU_time0();
-
cycles=E->control.p_iterations;
/* Solve for velocity and pressure, correct for bc's */
- residual_ddash=solve_Ahat_p_fhat(E,E->U,E->P,E->F,E->control.accuracy,&cycles);
+ solve_Ahat_p_fhat(E,E->U,E->P,E->F,E->control.accuracy,&cycles);
v_from_vector(E);
p_to_nodes(E,E->P,E->NP,E->mesh.levmax);
-/* */
-
- return;
+ return;
}
void solve_constrained_flow_iterative_pseudo_surf(E)
struct All_variables *E;
{
- double *D1;
- double *u;
- double *R,*Bp;
- double residual_ddash;
- double vmag;
- double global_vdot(),global_pdot();
-
float solve_Ahat_p_fhat();
void v_from_vector_pseudo_surf();
void p_to_nodes();
- int steps,cycles;
- int i,j,k,doff,vel_cycles_previous,vel_calls_previous;
+ int cycles;
- double time,CPU_time0();
-
- const int npno = E->lmesh.npno;
- const int gnpno = E->mesh.npno;
- const int nno = E->lmesh.nno;
- const int dims = E->mesh.nsd;
- const int neq = E->lmesh.neq;
- const int gneq = E->mesh.neq;
- const int addi_dof = additional_dof[dims];
-
- time=CPU_time0();
-
cycles=E->control.p_iterations;
/* Solve for velocity and pressure, correct for bc's */
- residual_ddash=solve_Ahat_p_fhat(E,E->U,E->P,E->F,E->control.accuracy,&cycles);
+ solve_Ahat_p_fhat(E,E->U,E->P,E->F,E->control.accuracy,&cycles);
v_from_vector_pseudo_surf(E);
p_to_nodes(E,E->P,E->NP,E->mesh.levmax);
-/* */
-
- return;
+ return;
}
-/* ========================================================================== */
+/* ========================================================================= */
-float solve_Ahat_p_fhat(E,V,P,FF,imp,steps_max)
-
- struct All_variables *E;
- double **V,**P,**FF;
- double imp;
- int *steps_max;
-
+float solve_Ahat_p_fhat(struct All_variables *E,
+ double **V, double **P, double **FF,
+ double imp, int *steps_max)
{
- int m,i,j,k,ii,count,convergent,valid,problems,lev,lev_low,npno,neq,steps;
- int gnpno,gneq;
+ int m, i, j, count, valid, lev, npno, neq;
+ int gnpno, gneq;
- double *r1[NCS],*F[NCS];
- double *r0[NCS],*r2[NCS],*z0[NCS],*z1[NCS],*s1[NCS],*s2[NCS],*Ah[NCS];
- double *shuffle[NCS];
- double alpha,delta,s2dotAhat,r0dotr0,r1dotz1;
- double residual, initial_residual, last_residual,v_res;
+ double *r1[NCS], *F[NCS];
+ double *r0[NCS], *r2[NCS], *z0[NCS], *z1[NCS], *s1[NCS], *s2[NCS];
+ double *shuffle[NCS];
+ double alpha, delta, s2dotAhat, r0dotr0, r1dotz1;
+ double residual, v_res;
- double global_vdot(),global_pdot();
- double *dvector();
+ double global_vdot(), global_pdot();
+ double *dvector();
- double time0,time,CPU_time0();
- float dpressure,dvelocity;
+ double time0, CPU_time0();
+ float dpressure, dvelocity;
- void assemble_div_u();
- void assemble_del2_u();
- void assemble_grad_p();
- void strip_bcs_from_residual();
- int solve_del2_u();
- void parallel_process_termination();
+ void assemble_div_u();
+ void assemble_del2_u();
+ void assemble_grad_p();
+ void strip_bcs_from_residual();
+ int solve_del2_u();
+ void parallel_process_termination();
- const int dims=E->mesh.nsd;
- const int n=loc_mat_size[E->mesh.nsd];
+ gnpno = E->mesh.npno;
+ gneq = E->mesh.neq;
+ npno = E->lmesh.npno;
+ neq = E->lmesh.neq;
- gnpno=E->mesh.npno;
- gneq=E->mesh.neq;
+ for (m=1; m<=E->sphere.caps_per_proc; m++) {
+ F[m] = (double *)malloc((neq+1)*sizeof(double));
- for (m=1;m<=E->sphere.caps_per_proc;m++) {
- npno=E->lmesh.npno;
- neq=E->lmesh.neq;
- F[m] = (double *)malloc((neq+1)*sizeof(double));
+ r0[m] = (double *)malloc((npno+1)*sizeof(double));
+ r1[m] = (double *)malloc((npno+1)*sizeof(double));
+ r2[m] = (double *)malloc((npno+1)*sizeof(double));
+ z0[m] = (double *)malloc((npno+1)*sizeof(double));
+ z1[m] = (double *)malloc((npno+1)*sizeof(double));
+ s1[m] = (double *)malloc((npno+1)*sizeof(double));
+ s2[m] = (double *)malloc((npno+1)*sizeof(double));
+ }
- r0[m] = (double *)malloc((npno+1)*sizeof(double));
- r1[m] = (double *)malloc((npno+1)*sizeof(double));
- r2[m] = (double *)malloc((npno+1)*sizeof(double));
- z0[m] = (double *)malloc((npno+1)*sizeof(double));
- z1[m] = (double *)malloc((npno+1)*sizeof(double));
- s1[m] = (double *)malloc((npno+1)*sizeof(double));
- s2[m] = (double *)malloc((npno+1)*sizeof(double));
- }
+ /* Copy the original force vector. FF shouldn't be modified. */
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ for(i=0;i<neq;i++)
+ F[m][i] = FF[m][i];
- /* Copy the original force vector. FF shouldn't be modified. */
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=0;i<neq;i++)
- F[m][i] = FF[m][i];
+ time0 = CPU_time0();
- problems=0;
- time0=time=CPU_time0();
+ /* calculate the initial velocity residual */
+ lev = E->mesh.levmax;
+ v_res = sqrt(global_vdot(E, F, F, lev)/gneq);
- /* calculate the velocity residual, note there are tricks involved here */
+ if (E->parallel.me==0) {
+ fprintf(E->fp, "initial residue of momentum equation F %.9e %d\n",
+ v_res, gneq);
+ fprintf(stderr, "initial residue of momentum equation F %.9e %d\n",
+ v_res, gneq);
+ }
- lev=E->mesh.levmax;
+ /* F = F - grad(P) - K*V */
+ assemble_grad_p(E, P, E->u1, lev);
+ for(m=1; m<=E->sphere.caps_per_proc; m++)
+ for(i=0; i<neq; i++)
+ F[m][i] = F[m][i] - E->u1[m][i];
- v_res=sqrt(global_vdot(E,F,F,lev)/gneq);
+ assemble_del2_u(E, V, E->u1, lev, 1);
+ for(m=1; m<=E->sphere.caps_per_proc; m++)
+ for(i=0; i<neq; i++)
+ F[m][i] = F[m][i] - E->u1[m][i];
+ strip_bcs_from_residual(E, F, lev);
- if (E->parallel.me==0) {
- fprintf(E->fp,"initial residue of momentum equation F %.9e %d\n",v_res,gneq);
- fprintf(stderr,"initial residue of momentum equation F %.9e %d\n",v_res,gneq);
- }
+ /* solve K*u1 = F for u1 */
+ valid=solve_del2_u(E, E->u1, F, imp*v_res, E->mesh.levmax);
+ strip_bcs_from_residual(E, E->u1, lev);
- assemble_grad_p(E,P,E->u1,lev);
+ /* V = V + u1 */
+ for(m=1; m<=E->sphere.caps_per_proc; m++)
+ for(i=0; i<neq; i++)
+ V[m][i] += E->u1[m][i];
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=0;i<neq;i++)
- F[m][i] -= E->u1[m][i];
- assemble_del2_u(E,V,E->u1,lev,1);
+ /* r1 = div(V) */
+ assemble_div_u(E, V, r1, lev);
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=0;i<neq;i++)
- F[m][i] -= E->u1[m][i];
+ /* incompressiblity residual = norm(r1) / norm(V) */
+ residual = sqrt(global_pdot(E, r1, r1, lev)/gnpno);
+ E->monitor.vdotv = sqrt(global_vdot(E, V, V, lev)/gneq);
+ E->monitor.incompressibility = residual / E->monitor.vdotv;
- strip_bcs_from_residual(E,F,lev);
+ count = 0;
- valid=solve_del2_u(E,E->u1,F,imp*v_res,E->mesh.levmax);
- strip_bcs_from_residual(E,E->u1,lev);
-
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=0;i<neq;i++)
- V[m][i] += E->u1[m][i];
-
- assemble_div_u(E,V,r1,lev);
-
- residual = initial_residual = sqrt(global_pdot(E,r1,r1,lev)/gnpno);
-
- E->monitor.vdotv = sqrt(global_vdot(E,V,V,lev)/gneq);
-
- E->monitor.incompressibility = residual/E->monitor.vdotv;
-
- count = 0;
- convergent=0;
-
- if (E->control.print_convergence && E->parallel.me==0) {
- fprintf(E->fp,"AhatP (%03d) after %g seconds with div/v=%.3e dv/v=%.3e"
- " and dp/p=%.3e for step %d\n",
- count, CPU_time0()-time0, E->monitor.incompressibility,
- 0.0, 0.0, E->monitor.solution_cycles);
- fflush(E->fp);
- fprintf(stderr,"AhatP (%03d) after %.3f seconds with v=%.3e div/v=%.3e dv/v=%.3e"
- " and dp/p=%.3e for step %d\n",
- count, CPU_time0()-time0, E->monitor.vdotv ,E->monitor.incompressibility,
- 0.0, 0.0, E->monitor.solution_cycles);
- }
-
- dpressure = 1.0;
- dvelocity = 1.0;
-
- while( (valid) && (count < *steps_max) &&
- (E->monitor.incompressibility >= E->control.tole_comp) &&
- (dpressure >= imp) && (dvelocity >= imp) ) {
-
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(j=1;j<=npno;j++)
- z1[m][j] = E->BPI[lev][m][j]*r1[m][j];
-
- r1dotz1 = global_pdot(E,r1,z1,lev);
-
- if ((count == 0))
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(j=1;j<=npno;j++)
- s2[m][j] = z1[m][j];
- else {
- r0dotr0=global_pdot(E,r0,z0,lev);
- assert(r0dotr0 != 0.0 /* Division by zero in head of incompressibility iteration */);
- delta = r1dotz1/r0dotr0;
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(j=1;j<=npno;j++)
- s2[m][j] = z1[m][j] + delta * s1[m][j];
+ if (E->control.print_convergence && E->parallel.me==0) {
+ fprintf(E->fp, "AhatP (%03d) after %g seconds with div/v=%.3e "
+ "for step %d\n", count, CPU_time0()-time0,
+ E->monitor.incompressibility, E->monitor.solution_cycles);
+ fprintf(stderr, "AhatP (%03d) after %g seconds with div/v=%.3e "
+ "for step %d\n", count, CPU_time0()-time0,
+ E->monitor.incompressibility, E->monitor.solution_cycles);
}
- assemble_grad_p(E,s2,F,lev);
- valid=solve_del2_u(E,E->u1,F,imp*v_res,lev);
- strip_bcs_from_residual(E,E->u1,lev);
+ /* pressure and velocity corrections */
+ dpressure = 1.0;
+ dvelocity = 1.0;
- assemble_div_u(E,E->u1,F,lev);
+ while( (valid) && (count < *steps_max) &&
+ (E->monitor.incompressibility >= E->control.tole_comp) &&
+ (dpressure >= imp) && (dvelocity >= imp) ) {
- s2dotAhat=global_pdot(E,s2,F,lev);
- if (valid)
- /* alpha defined this way is the same as R&W */
- alpha = r1dotz1/s2dotAhat;
- else
- alpha = 0.0;
+ /* preconditioner B, z1 = B*r1 */
+ for(m=1; m<=E->sphere.caps_per_proc; m++)
+ for(j=1; j<=npno; j++)
+ z1[m][j] = E->BPI[lev][m][j] * r1[m][j];
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(j=1;j<=npno;j++) {
- r2[m][j] = r1[m][j] - alpha * F[m][j];
- P[m][j] += alpha * s2[m][j];
- }
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(j=0;j<neq;j++)
- V[m][j] -= alpha * E->u1[m][j];
+ /* r1dotz1 = <r1, z1> */
+ r1dotz1 = global_pdot(E, r1, z1, lev);
- assemble_div_u(E,V,F,lev);
- E->monitor.vdotv = global_vdot(E,V,V,E->mesh.levmax);
- E->monitor.incompressibility = sqrt((gneq/gnpno)*(1.0e-32+global_pdot(E,F,F,lev)/(1.0e-32+E->monitor.vdotv)));
- dpressure = alpha * sqrt(global_pdot(E,s2,s2,lev)
- / (1.0e-32 + global_pdot(E,P,P,lev)));
- dvelocity = alpha * sqrt(global_vdot(E,E->u1,E->u1,lev)
- / (1.0e-32 + E->monitor.vdotv));
+ if ((count == 0))
+ for (m=1; m<=E->sphere.caps_per_proc; m++)
+ for(j=1; j<=npno; j++)
+ s2[m][j] = z1[m][j];
+ else {
+ /* s2 = z1 + s1 * <r1,z1>/<r0,z0> */
+ r0dotr0 = global_pdot(E, r0, z0, lev);
+ assert(r0dotr0 != 0.0 /* Division by zero in head of incompressibility iteration */);
+ delta = r1dotz1 / r0dotr0;
+ for(m=1; m<=E->sphere.caps_per_proc; m++)
+ for(j=1; j<=npno; j++)
+ s2[m][j] = z1[m][j] + delta * s1[m][j];
+ }
- count++;
- if (E->control.print_convergence && E->parallel.me==0) {
- fprintf(E->fp, "AhatP (%03d) after %g seconds with div/v=%.3e dv/v=%.3e"
- " and dp/p=%.3e for step %d\n",
- count, CPU_time0()-time0, E->monitor.incompressibility,
- dvelocity, dpressure, E->monitor.solution_cycles);
- fprintf(stderr, "AhatP (%03d) after %.3f seconds with v=%.3e div/v=%.3e dv/v=%.3e"
- " and dp/p=%.3e for step %d\n",
- count, CPU_time0()-time0, E->monitor.vdotv,E->monitor.incompressibility,
- dvelocity, dpressure, E->monitor.solution_cycles);
- }
+ /* solve K*u1 = grad(s2) for u1 */
+ assemble_grad_p(E, s2, F, lev);
+ valid = solve_del2_u(E, E->u1, F, imp*v_res, lev);
+ strip_bcs_from_residual(E, E->u1, lev);
- for (m=1;m<=E->sphere.caps_per_proc;m++) {
- shuffle[m]=s1[m];s1[m]=s2[m];s2[m]=shuffle[m];
- shuffle[m]=r0[m];r0[m]=r1[m];r1[m]=r2[m];r2[m]=shuffle[m];
- shuffle[m]=z0[m];z0[m]=z1[m];z1[m]=shuffle[m];
- }
- } /* end loop for conjugate gradient */
+ /* alpha = <r1, z1> / <s2, div(u1)> */
+ assemble_div_u(E, E->u1, F, lev);
+ s2dotAhat = global_pdot(E, s2, F, lev);
- if(problems) {
- fprintf(E->fp,"Convergence of velocity solver may affect continuity\n");
- fprintf(E->fp,"Consider running with the `see_convergence=on' option\n");
- fprintf(E->fp,"To evaluate the performance of the current relaxation parameters\n");
- fflush(E->fp);
- }
+ if(valid)
+ /* alpha defined this way is the same as R&W */
+ alpha = r1dotz1 / s2dotAhat;
+ else
+ alpha = 0.0;
- for (m=1;m<=E->sphere.caps_per_proc;m++) {
- free((void *) F[m]);
- free((void *) r0[m]);
- free((void *) r1[m]);
- free((void *) r2[m]);
- free((void *) z0[m]);
- free((void *) z1[m]);
- free((void *) s1[m]);
- free((void *) s2[m]);
- }
- *steps_max=count;
+ /* r2 = r1 - alpha * div(u1) */
+ /* P = P + alpha * s2 */
+ /* V = V - alpha * u1 */
+ for(m=1; m<=E->sphere.caps_per_proc; m++)
+ for(j=1; j<=npno; j++) {
+ r2[m][j] = r1[m][j] - alpha * F[m][j];
+ P[m][j] += alpha * s2[m][j];
+ }
- return(residual);
-}
+ for(m=1; m<=E->sphere.caps_per_proc; m++)
+ for(j=0; j<neq; j++)
+ V[m][j] -= alpha * E->u1[m][j];
-/* ========================================================================== */
+ /* compute velocity and incompressibility residual */
+ assemble_div_u(E, V, F, lev);
+ E->monitor.vdotv = global_vdot(E, V, V, E->mesh.levmax);
+ E->monitor.incompressibility = sqrt((gneq/gnpno)
+ *(1.0e-32
+ + global_pdot(E, F, F, lev)
+ / (1.0e-32+E->monitor.vdotv)));
+ /* compute velocity and pressure corrections */
+ dpressure = alpha * sqrt(global_pdot(E, s2, s2, lev)
+ / (1.0e-32 + global_pdot(E, P, P, lev)));
+ dvelocity = alpha * sqrt(global_vdot(E, E->u1, E->u1, lev)
+ / (1.0e-32 + E->monitor.vdotv));
-void v_from_vector(E)
- struct All_variables *E;
-{
- int m,node;
- const int nno = E->lmesh.nno;
+ count++;
- for (m=1;m<=E->sphere.caps_per_proc;m++) {
- for(node=1;node<=nno;node++) {
- E->sphere.cap[m].V[1][node] = E->U[m][E->id[m][node].doff[1]];
- E->sphere.cap[m].V[2][node] = E->U[m][E->id[m][node].doff[2]];
- E->sphere.cap[m].V[3][node] = E->U[m][E->id[m][node].doff[3]];
- if (E->node[m][node] & VBX)
- E->sphere.cap[m].V[1][node] = E->sphere.cap[m].VB[1][node];
- if (E->node[m][node] & VBY)
- E->sphere.cap[m].V[2][node] = E->sphere.cap[m].VB[2][node];
- if (E->node[m][node] & VBZ)
- E->sphere.cap[m].V[3][node] = E->sphere.cap[m].VB[3][node];
- }
- }
+ if(E->control.print_convergence && E->parallel.me==0) {
+ fprintf(E->fp, "AhatP (%03d) after %g seconds with div/v=%.3e "
+ "dv/v=%.3e and dp/p=%.3e for step %d\n",
+ count, CPU_time0()-time0, E->monitor.incompressibility,
+ dvelocity, dpressure, E->monitor.solution_cycles);
+ fprintf(stderr, "AhatP (%03d) after %g seconds with div/v=%.3e "
+ "dv/v=%.3e and dp/p=%.3e for step %d\n",
+ count, CPU_time0()-time0, E->monitor.incompressibility,
+ dvelocity, dpressure, E->monitor.solution_cycles);
+ }
- return;
-}
-void v_from_vector_pseudo_surf(E)
- struct All_variables *E;
-{
- int m,node;
+ /* swap array pointers */
+ for(m=1; m<=E->sphere.caps_per_proc; m++) {
+ shuffle[m] = s1[m];
+ s1[m] = s2[m];
+ s2[m] = shuffle[m];
- const int nno = E->lmesh.nno;
- const int dofs = E->mesh.dof;
- double sum_V = 0.0, sum_dV = 0.0, rel_error = 0.0, global_max_error = 0.0;
- double tol_error = 1.0e-03;
+ shuffle[m] = r0[m];
+ r0[m] = r1[m];
+ r1[m] = r2[m];
+ r2[m] = shuffle[m];
- for (m=1;m<=E->sphere.caps_per_proc;m++) {
- for(node=1;node<=nno;node++) {
- E->sphere.cap[m].Vprev[1][node] = E->sphere.cap[m].V[1][node];
- E->sphere.cap[m].Vprev[2][node] = E->sphere.cap[m].V[2][node];
- E->sphere.cap[m].Vprev[3][node] = E->sphere.cap[m].V[3][node];
+ shuffle[m] = z0[m];
+ z0[m] = z1[m];
+ z1[m] = shuffle[m];
+ }
- E->sphere.cap[m].V[1][node] = E->U[m][E->id[m][node].doff[1]];
- E->sphere.cap[m].V[2][node] = E->U[m][E->id[m][node].doff[2]];
- E->sphere.cap[m].V[3][node] = E->U[m][E->id[m][node].doff[3]];
- if (E->node[m][node] & VBX)
- E->sphere.cap[m].V[1][node] = E->sphere.cap[m].VB[1][node];
- if (E->node[m][node] & VBY)
- E->sphere.cap[m].V[2][node] = E->sphere.cap[m].VB[2][node];
- if (E->node[m][node] & VBZ)
- E->sphere.cap[m].V[3][node] = E->sphere.cap[m].VB[3][node];
+ } /* end loop for conjugate gradient */
- sum_dV += (E->sphere.cap[m].V[1][node] - E->sphere.cap[m].Vprev[1][node])*(E->sphere.cap[m].V[1][node] - E->sphere.cap[m].Vprev[1][node])
- + (E->sphere.cap[m].V[2][node] - E->sphere.cap[m].Vprev[2][node])*(E->sphere.cap[m].V[2][node] - E->sphere.cap[m].Vprev[2][node])
- + (E->sphere.cap[m].V[3][node] - E->sphere.cap[m].Vprev[3][node])*(E->sphere.cap[m].V[3][node] - E->sphere.cap[m].Vprev[3][node]);
- sum_V += E->sphere.cap[m].V[1][node]*E->sphere.cap[m].V[1][node]
- + E->sphere.cap[m].V[2][node]*E->sphere.cap[m].V[2][node]
- + E->sphere.cap[m].V[3][node]*E->sphere.cap[m].V[3][node];
- }
- rel_error = sqrt(sum_dV)/sqrt(sum_V);
- MPI_Allreduce(&rel_error,&global_max_error,1,MPI_DOUBLE,MPI_MAX,E->parallel.world);
- if(global_max_error <= tol_error) E->monitor.stop_topo_loop = 1;
- if(E->parallel.me==0)
- fprintf(stderr,"global_max_error=%e stop_topo_loop=%d\n",global_max_error,E->monitor.stop_topo_loop);
+ for(m=1; m<=E->sphere.caps_per_proc; m++) {
+ free((void *) r0[m]);
+ free((void *) r1[m]);
+ free((void *) r2[m]);
+ free((void *) z0[m]);
+ free((void *) z1[m]);
+ free((void *) s1[m]);
+ free((void *) s2[m]);
+ }
- }
+ *steps_max=count;
- return;
+ return(residual);
}
-void velo_from_element(E,VV,m,el,sphere_key)
- struct All_variables *E;
- float VV[4][9];
- int el,m,sphere_key;
- {
+/* ========================================================================== */
- int a, node;
- double sint, cost, sinf, cosf;
- const int dims=E->mesh.nsd;
- const int ends=enodes[E->mesh.nsd];
- const int nno=E->lmesh.nno;
- const int lev=E->mesh.levmax;
- if (sphere_key)
- for(a=1;a<=ends;a++) {
- node = E->ien[m][el].node[a];
- VV[1][a] = E->sphere.cap[m].V[1][node];
- VV[2][a] = E->sphere.cap[m].V[2][node];
- VV[3][a] = E->sphere.cap[m].V[3][node];
- }
- else {
- for(a=1;a<=ends;a++) {
- node = E->ien[m][el].node[a];
- sint = E->SinCos[lev][m][0][node];
- sinf = E->SinCos[lev][m][1][node];
- cost = E->SinCos[lev][m][2][node];
- cosf = E->SinCos[lev][m][3][node];
- VV[1][a] = E->sphere.cap[m].V[1][node]*cost*cosf
- - E->sphere.cap[m].V[2][node]*sinf
- + E->sphere.cap[m].V[3][node]*sint*cosf;
- VV[2][a] = E->sphere.cap[m].V[1][node]*cost*sinf
- + E->sphere.cap[m].V[2][node]*cosf
- + E->sphere.cap[m].V[3][node]*sint*sinf;
- VV[3][a] = -E->sphere.cap[m].V[1][node]*sint
- + E->sphere.cap[m].V[3][node]*cost;
- }
- }
- return;
- }
-
-
-void velo_from_element_d(E,VV,m,el,sphere_key)
- struct All_variables *E;
- double VV[4][9];
- int el,m,sphere_key;
-{
-
- int a, node;
- double sint, cost, sinf, cosf;
- const int dims=E->mesh.nsd;
- const int ends=enodes[E->mesh.nsd];
- const int nno=E->lmesh.nno;
- const int lev=E->mesh.levmax;
-
- if (sphere_key)
- for(a=1;a<=ends;a++) {
- node = E->ien[m][el].node[a];
- VV[1][a] = E->sphere.cap[m].V[1][node];
- VV[2][a] = E->sphere.cap[m].V[2][node];
- VV[3][a] = E->sphere.cap[m].V[3][node];
- }
- else {
- for(a=1;a<=ends;a++) {
- node = E->ien[m][el].node[a];
-
- sint = E->SinCos[lev][m][0][node];
- sinf = E->SinCos[lev][m][1][node];
- cost = E->SinCos[lev][m][2][node];
- cosf = E->SinCos[lev][m][3][node];
-
- VV[1][a] = E->sphere.cap[m].V[1][node]*cost*cosf
- - E->sphere.cap[m].V[2][node]*sinf
- + E->sphere.cap[m].V[3][node]*sint*cosf;
- VV[2][a] = E->sphere.cap[m].V[1][node]*cost*sinf
- + E->sphere.cap[m].V[2][node]*cosf
- + E->sphere.cap[m].V[3][node]*sint*sinf;
- VV[3][a] = -E->sphere.cap[m].V[1][node]*sint
- + E->sphere.cap[m].V[3][node]*cost;
- }
- }
- return;
-}
More information about the cig-commits
mailing list