Revision 408e100cc1e612db73f405c28463f7606dee3090 authored by Eh Tan on 10 August 2007, 23:56:37 UTC, committed by Eh Tan on 10 August 2007, 23:56:37 UTC
1 parent d603970
misc.cc
// -*- C++ -*-
//
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
//<LicenseText>
//
// CitcomS.py by Eh Tan, Eun-seo Choi, and Pururav Thoutireddy.
// Copyright (C) 2002-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>
//
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
#include "config.h"
#include <Python.h>
//#include <iostream>
//#include <fstream>
#include "mpi.h"
#include "mpi/pympi.h"
#include "misc.h"
extern "C" {
#include "global_defs.h"
#include "element_definitions.h"
#include "citcom_init.h"
}
void commonE(All_variables*);
// copyright
char PyCitcomSExchanger_copyright__doc__[] = "";
char PyCitcomSExchanger_copyright__name__[] = "copyright";
static char PyCitcomSExchanger_copyright_note[] =
"CitcomS Exchanger python module: Copyright (c) 1998-2003 California Institute of Technology";
PyObject * PyCitcomSExchanger_copyright(PyObject *, PyObject *)
{
return Py_BuildValue("s", PyCitcomSExchanger_copyright_note);
}
// return (All_variables* E)
char PyCitcomSExchanger_FinereturnE__doc__[] = "";
char PyCitcomSExchanger_FinereturnE__name__[] = "FinereturnE";
PyObject * PyCitcomSExchanger_FinereturnE(PyObject *, PyObject *args)
{
PyObject *Obj;
if (!PyArg_ParseTuple(args, "O:FinereturnE", &Obj))
return NULL;
PyMPICommObject* comm = (PyMPICommObject *)Obj;
MPI_Comm world = comm->comm;
All_variables *E = citcom_init(&world);
E->lmesh.nox = 4;
E->lmesh.noy = 4;
E->lmesh.noz = 3;
E->control.theta_max = 2.1;
E->control.theta_min = 0.9;
E->control.fi_max = 2.1;
E->control.fi_min = 0.9;
E->sphere.ro = 1.9;
E->sphere.ri = 0.9;
commonE(E);
PyObject *cobj = PyCObject_FromVoidPtr(E, NULL);
return Py_BuildValue("O", cobj);
}
char PyCitcomSExchanger_CoarsereturnE__doc__[] = "";
char PyCitcomSExchanger_CoarsereturnE__name__[] = "CoarsereturnE";
PyObject * PyCitcomSExchanger_CoarsereturnE(PyObject *, PyObject *args)
{
PyObject *Obj;
if (!PyArg_ParseTuple(args, "O:CoarsereturnE", &Obj))
return NULL;
PyMPICommObject* comm = (PyMPICommObject *)Obj;
MPI_Comm world = comm->comm;
All_variables *E = citcom_init(&world);
E->lmesh.nox = 4;
E->lmesh.noy = 4;
E->lmesh.noz = 3;
E->control.theta_max = 3.0;
E->control.theta_min = 0.0;
E->control.fi_max = 3.0;
E->control.fi_min = 0.0;
E->sphere.ro = 2.0;
E->sphere.ri = 0.0;
commonE(E);
// std::ofstream cfile("coarse.dat");
for(int m=1;m<=E->sphere.caps_per_proc;m++)
for(int k=1;k<=E->lmesh.noy;k++)
for(int j=1;j<=E->lmesh.nox;j++)
for(int i=1;i<=E->lmesh.noz;i++) {
int node = i + (j-1)*E->lmesh.noz
+ (k-1)*E->lmesh.noz*E->lmesh.nox;
E->T[m][node] = E->sx[m][1][node]
+ E->sx[m][2][node]
+ E->sx[m][3][node];
E->sphere.cap[m].V[1][node] = E->T[m][node];
E->sphere.cap[m].V[2][node] = 2.0*E->T[m][node];
E->sphere.cap[m].V[3][node] = 3.0*E->T[m][node];
// cfile << "in CoarsereturnE (T, v1,v2,v3): "
// << node << " "
// << E->sx[m][1][node] << " "
// << E->sx[m][2][node] << " "
// << E->sx[m][3][node] << " "
// << E->T[m][node] << " "
// << E->sphere.cap[m].V[1][node] << " "
// << E->sphere.cap[m].V[2][node] << " "
// << E->sphere.cap[m].V[3][node] << " "
// << std::endl;
}
// cfile.close();
PyObject *cobj = PyCObject_FromVoidPtr(E, NULL);
return Py_BuildValue("O", cobj);
}
void commonE(All_variables *E)
{
E->control.accuracy = 1e-6;
E->control.tole_comp = 1e-7;
E->parallel.nprocxy = 1;
E->parallel.nprocx = E->parallel.nproc;
E->parallel.nprocy = 1;
E->parallel.nprocz = 1;
E->parallel.me_loc[1] = E->parallel.me;
E->parallel.me_loc[2] = 0;
E->parallel.me_loc[3] = 0;
E->sphere.caps_per_proc = 1;
E->mesh.levmax = 1;
E->mesh.levmin = 1;
E->mesh.dof = 3;
E->lmesh.elx = E->lmesh.nox - 1;
E->lmesh.elz = E->lmesh.noz - 1;
E->lmesh.ely = E->lmesh.noy - 1;
E->lmesh.nno = E->lmesh.noz * E->lmesh.nox * E->lmesh.noy;
E->lmesh.nel = E->lmesh.ely * E->lmesh.elx * E->lmesh.elz;
E->lmesh.npno = E->lmesh.nel;
int noz = E->lmesh.noz;
int noy = E->lmesh.noy;
int nox = E->lmesh.nox;
E->lmesh.ELX[E->mesh.levmax] = nox-1;
E->lmesh.ELY[E->mesh.levmax] = noy-1;
E->lmesh.ELZ[E->mesh.levmax] = noz-1;
E->lmesh.NOZ[E->mesh.levmax] = noz;
E->lmesh.NOY[E->mesh.levmax] = noy;
E->lmesh.NOX[E->mesh.levmax] = nox;
E->lmesh.NNO[E->mesh.levmax] = nox * noz * noy;
E->lmesh.NEL[E->mesh.levmax] = (nox-1) * (noz-1) * (noy-1);
E->mesh.elx = E->lmesh.elx * E->parallel.nprocx;
E->mesh.elz = E->lmesh.elz * E->parallel.nprocz;
E->mesh.ely = E->lmesh.ely * E->parallel.nprocy;
E->mesh.nox = E->mesh.elx + 1;
E->mesh.noz = E->mesh.elz + 1;
E->mesh.noy = E->mesh.ely + 1;
E->mesh.nno = E->mesh.nox * E->mesh.noy * E->mesh.noz;
E->mesh.nel = E->mesh.ely * E->mesh.elx * E->mesh.elz;
for (int j=1;j<=E->sphere.caps_per_proc;j++) {
E->sphere.capid[j] = 1;
}
for (int lev=E->mesh.levmax;lev>=E->mesh.levmin;lev--) {
for (int j=1;j<=E->sphere.caps_per_proc;j++) {
E->IEN[lev][j] = new IEN [E->lmesh.nel+1];
E->ECO[lev][j] = (struct COORD *) malloc((E->lmesh.nno+2)*sizeof(struct COORD));
}
}
for (int j=1;j<=E->sphere.caps_per_proc;j++) {
E->ien[j] = E->IEN[E->mesh.levmax][j];
E->eco[j] = E->ECO[E->mesh.levmax][j];
}
for (int lev=E->mesh.levmax;lev>=E->mesh.levmin;lev--)
for (int j=1;j<=E->sphere.caps_per_proc;j++)
for(int n=1; n<=E->lmesh.nno; n++) {
E->ECO[lev][j][n].area = 1.0;
}
for(int m=1;m<=E->sphere.caps_per_proc;m++) {
E->sphere.cap[m].TB[1] = new float[E->lmesh.nno+1];
E->sphere.cap[m].TB[2] = new float[E->lmesh.nno+1];
E->sphere.cap[m].TB[3] = new float[E->lmesh.nno+1];
E->sphere.cap[m].VB[1] = new float[E->lmesh.nno+1];
E->sphere.cap[m].VB[2] = new float[E->lmesh.nno+1];
E->sphere.cap[m].VB[3] = new float[E->lmesh.nno+1];
E->node[m] = new unsigned int[E->lmesh.nno+1];
for(int n=1; n<=E->lmesh.nno; n++) {
E->sphere.cap[m].TB[1][n] = 0;
E->sphere.cap[m].TB[2][n] = 0;
E->sphere.cap[m].TB[3][n] = 0;
E->sphere.cap[m].VB[1][n] = 0;
E->sphere.cap[m].VB[2][n] = 0;
E->sphere.cap[m].VB[3][n] = 0;
E->node[m][n] = 0;
}
}
for (int lev=E->mesh.levmax;lev>=E->mesh.levmin;lev--) {
for (int j=1;j<=E->sphere.caps_per_proc;j++) {
int elx = E->lmesh.ELX[lev];
int elz = E->lmesh.ELZ[lev];
int ely = E->lmesh.ELY[lev];
int nox = E->lmesh.NOX[lev];
int noz = E->lmesh.NOZ[lev];
for(int r=1;r<=ely;r++)
for(int q=1;q<=elx;q++)
for(int p=1;p<=elz;p++) {
int element = (r-1)*elx*elz + (q-1)*elz + p;
int start = (r-1)*noz*nox + (q-1)*noz + p;
for(int rr=1;rr<=8;rr++) {
E->IEN[lev][j][element].node[rr]= start
+ offset[rr].vector[0]
+ offset[rr].vector[1]*noz
+ offset[rr].vector[2]*noz*nox;
// std::cout << " el = " << element
// << " lnode = " << rr
// << " ien = "
// << E->IEN[lev][j][element].node[rr]
// << std::endl;
}
}
} /* end for cap j */
} /* end loop for lev */
const int n = E->lmesh.nno;
for(int m=1;m<=E->sphere.caps_per_proc;m++) {
for(int i=1; i<=E->mesh.dof; i++) {
// Don't forget to delete these later
E->sx[m][i] = new double [n+1];
E->sphere.cap[m].V[i] = new float [n+1];
}
E->T[m] = new double [n+1];
}
for(int m=1;m<=E->sphere.caps_per_proc;m++)
for(int k=1;k<=E->lmesh.noy;k++)
for(int j=1;j<=E->lmesh.nox;j++)
for(int i=1;i<=E->lmesh.noz;i++) {
int node = i + (j-1)*E->lmesh.noz
+ (k-1)*E->lmesh.noz*E->lmesh.nox;
E->sx[m][1][node] =
(E->control.theta_max - E->control.theta_min)
/ E->mesh.elx * (j-1)
+ (E->control.theta_max - E->control.theta_min)
* E->parallel.me_loc[1] / E->parallel.nprocx
+ E->control.theta_min;
E->sx[m][2][node] =
(E->control.fi_max - E->control.fi_min)
/ E->mesh.ely * (k-1)
+ (E->control.fi_max - E->control.fi_min)
* E->parallel.me_loc[2] / E->parallel.nprocy
+ E->control.fi_min;
E->sx[m][3][node] =
(E->sphere.ro - E->sphere.ri)
/ E->mesh.elz * (i-1)
+ (E->sphere.ro - E->sphere.ri)
* E->parallel.me_loc[3] / E->parallel.nprocz
+ E->sphere.ri;
// std::cout << node << " "
// << E->sx[m][1][node] << " "
// << E->sx[m][2][node] << " "
// << E->sx[m][3][node] << " "
// << std::endl;
}
return;
}
// version
// $Id$
// End of file
Computing file changes ...