// -*- C++ -*-
//
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
//
//
// 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
//
//
//
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
#include
#include
#include
#include
#include "global_bbox.h"
#include "global_defs.h"
#include "journal/diagnostics.h"
#include "Exchanger/BoundedBox.h"
#include "Exchanger/BoundedMesh.h"
#include "CitcomInterpolator.h"
using Exchanger::Array2D;
using Exchanger::BoundedBox;
using Exchanger::BoundedMesh;
using Exchanger::DIM;
using Exchanger::NODES_PER_ELEMENT;
using Exchanger::STRESS_DIM;
CitcomInterpolator::CitcomInterpolator(const BoundedMesh& boundedMesh,
Array2D& meshNode,
const All_variables* e) :
Exchanger::Interpolator(),
E(e)
{
journal::debug_t debug("CitcomS-Exchanger");
debug << journal::at(__HERE__) << journal::endl;
init(boundedMesh, meshNode);
selfTest(boundedMesh, meshNode);
elem_.print("CitcomS-CitcomInterpolator-elem");
shape_.print("CitcomS-CitcomInterpolator-shape");
}
CitcomInterpolator::~CitcomInterpolator()
{}
void CitcomInterpolator::interpolatePressure(Array2D& P)
{
P.assign(size(), 0);
const int mm = 1;
for(int i=0; iien[mm][n1].node[k+1];
P[0][i] += shape_[k][i] * E->P[mm][node];
}
}
}
void CitcomInterpolator::interpolateStress(Array2D& S)
{
S.assign(size(), 0);
const int mm = 1;
for(int i=0; iien[mm][n1].node[k+1] - 1;
for(int d=0; dgstress[mm][node*STRESS_DIM+d+1];
}
}
}
void CitcomInterpolator::interpolateTemperature(Array2D& T)
{
T.assign(size(), 0);
const int mm = 1;
for(int i=0; iien[mm][n1].node[k+1];
T[0][i] += shape_[k][i] * E->T[mm][node];
}
}
}
void CitcomInterpolator::interpolateVelocity(Array2D& V)
{
V.assign(size(), 0);
const int mm = 1;
for(int i=0; iien[mm][n1].node[k+1];
for(int d=0; dsphere.cap[mm].V[d+1][node];
}
}
}
// private functions
void CitcomInterpolator::init(const BoundedMesh& boundedMesh,
Array2D& meshNode)
{
journal::debug_t debug("CitcomS-Exchanger");
debug << journal::at(__HERE__) << journal::endl;
const int mm = 1;
elem_.reserve(boundedMesh.size());
shape_.reserve(boundedMesh.size());
Array2D etaAxes; // axes of eta coord.
Array2D inv_length_sq; // reciprocal of (length of etaAxes)^2
computeElementGeometry(etaAxes, inv_length_sq);
// get remote BoundedBox
BoundedBox remoteBBox(boundedMesh.tightBBox());
// get local BoundedBox
BoundedBox bbox(DIM);
if(E->parallel.nprocxy == 12) {
// for CitcomS Full
fullGlobalBoundedBox(bbox, E);
}
else {
// for CitcomS Regional
regionalGlobalBoundedBox(bbox, E);
}
// z is the range of depth in current processor
std::vector z(E->lmesh.noz);
for(size_t i=0; isx[mm][DIM][i+1];
for(int n=0; n x(DIM);
for(int d=0; d 0)
x[DIM-1] = bbox[1][DIM-1];
}
// skip if x is not inside bbox
if(!isInside(x, bbox)) continue;
#if 1
// skip if x is not in the range of depth
if(x.back() < z.front() || x.back() > z.back()) continue;
int elz = bisectInsertPoint(x.back(), z);
// Since the mesh of CitcomS is structural and regular
// we only need to loop over elements in a constant depth
for(int el=elz; ellmesh.nel; el+=E->lmesh.elz) {
#else
// If the mesh is not regular, loop over all elements
for(int el=0; ellmesh.nel; ++el) {
#endif
std::vector elmShape(NODES_PER_ELEMENT);
double accuracy = E->control.accuracy * E->control.accuracy
* E->eco[mm][el+1].area;
bool found = elementInverseMapping(elmShape, x,
etaAxes, inv_length_sq,
el, accuracy);
if(found) {
meshNode.push_back(n);
elem_.push_back(el+1);
shape_.push_back(elmShape);
break;
}
}
}
}
void CitcomInterpolator::computeElementGeometry(Array2D& etaAxes,
Array2D& inv_length_sq) const
{
etaAxes.resize(E->lmesh.nel);
inv_length_sq.resize(E->lmesh.nel);
const int mm = 1;
const int surfnodes = 4; // # of nodes on element's face
// node 1, 2, 5, 6 are in the face that is on positive x axis
// node 2, 3, 6, 7 are in the face that is on positive y axis
// node 4, 5, 6, 7 are in the face that is on positive z axis
// node 0, 3, 4, 7 are in the face that is on negative x axis
// node 0, 1, 4, 5 are in the face that is on negative y axis
// node 0, 1, 2, 3 are in the face that is on negative z axis
// see comment of getShapeFunction() for the ordering of nodes
const int high[] = {1, 2, 5, 6,
2, 3, 6, 7,
4, 5, 6, 7};
const int low[] = {0, 3, 4, 7,
0, 1, 4, 5,
0, 1, 2, 3};
std::vector node(NODES_PER_ELEMENT);
for(int element=0; elementlmesh.nel; ++element) {
for(int k=0; kien[mm][element+1].node[k+1];
for(int n=0; nsx[mm][d+1][node[high[k+i]]];
lowmean += E->sx[mm][d+1][node[low[k+i]]];
}
etaAxes[n*DIM+d][element] = 0.5 * (highmean - lowmean)
/ surfnodes;
}
}
}
for(int element=0; elementlmesh.nel; ++element)
for(int n=0; n& v) const
{
int low = 0;
int high = v.size();
int insert_point = (low + high) / 2;
while(low < high) {
if(x < v[insert_point])
high = insert_point;
else if(x > v[insert_point+1])
low = insert_point;
else
break;
insert_point = (low + high) / 2;
}
return insert_point;
}
bool CitcomInterpolator::elementInverseMapping(std::vector& elmShape,
const std::vector& x,
const Array2D& etaAxes,
const Array2D& inv_length_sq,
int element,
double accuracy)
{
const int mm = 1;
bool found = false;
int count = 0;
std::vector eta(DIM); // initial eta = (0,0,0)
while (1) {
getShapeFunction(elmShape, eta);
std::vector xx(DIM);
for(int k=0; kien[mm][element+1].node[k+1];
for(int d=0; dsx[mm][d+1][node] * elmShape[k];
}
std::vector dx(DIM);
double distancesq = 0;
for(int d=0; d deta(DIM);
for(int d=0; d 1.5);
// iterate at least twice
found = (distancesq < accuracy) && (count > 0);
if (outside || found || (count > 100))
break;
++count;
/* Only need to iterate if this is marginal. If eta > distortion of
an individual element then almost certainly x is in a
different element ... or the mesh is terrible ! */
}
return found;
}
void CitcomInterpolator::getShapeFunction(std::vector& shape,
const std::vector& eta) const
{
// the ordering of nodes in an element
// node #: eta coordinate
// node 0: (-1, -1, -1)
// node 1: ( 1, -1, -1)
// node 2: ( 1, 1, -1)
// node 3: (-1, 1, -1)
// node 4: (-1, -1, 1)
// node 5: ( 1, -1, 1)
// node 6: ( 1, 1, 1)
// node 7: (-1, 1, 1)
shape[0] = 0.125 * (1.0-eta[0]) * (1.0-eta[1]) * (1.0-eta[2]);
shape[1] = 0.125 * (1.0+eta[0]) * (1.0-eta[1]) * (1.0-eta[2]);
shape[2] = 0.125 * (1.0+eta[0]) * (1.0+eta[1]) * (1.0-eta[2]);
shape[3] = 0.125 * (1.0-eta[0]) * (1.0+eta[1]) * (1.0-eta[2]);
shape[4] = 0.125 * (1.0-eta[0]) * (1.0-eta[1]) * (1.0+eta[2]);
shape[5] = 0.125 * (1.0+eta[0]) * (1.0-eta[1]) * (1.0+eta[2]);
shape[6] = 0.125 * (1.0+eta[0]) * (1.0+eta[1]) * (1.0+eta[2]);
shape[7] = 0.125 * (1.0-eta[0]) * (1.0+eta[1]) * (1.0+eta[2]);
}
void CitcomInterpolator::selfTest(const BoundedMesh& boundedMesh,
const Array2D& meshNode) const
{
journal::debug_t debug("CitcomS-Exchanger");
debug << journal::at(__HERE__) << journal::endl;
for(int i=0; i xt(DIM);
for(int j=0; j xi(DIM);
int n1 = elem_[0][i];
for(int j=0; jien[1][n1].node[j+1];
for(int k=0; ksx[1][k+1][node] * shape_[j][i];
}
// if xi and xt are not coincide, the interpolation is wrong
double norm = 0.0;
for(int k=0; k 1.e-10) {
double tshape = 0.0;
for(int j=0; j