Commit cfb923f1 authored by thomas's avatar thomas
Browse files

dg: GBR testcase update

git-svn-id: https://geuz.org/svn/gmsh/trunk@13555 59c8cabf-10d6-4a28-b43f-3936b85956c4
parent 6b95e9a8
......@@ -355,9 +355,11 @@ else:
# ***SEED AT A POINT*** (NB: Look at method used to print results if changing seeding method)
#cartCoords = [0]*3
#geoCoordToCartesian(-19.9400,149.1100,cartCoords)
#geoCoordTo3DCartesian(-19.9400,149.1100,cartCoords)
#projectPoint(cartCoords)
#particles = particleArray(groups, 1000, cartCoords[0], cartCoords[1], 0)
#particles = particleArray(groups, 0, cartCoords[0], cartCoords[1], 0, 0)
#for location in range(0,10):
#particles.addParticles(groups, 1000, cartCoords[0], cartCoords[1], 0, location)
#particles.printPositions(LP_OutputDir+'/seededParticles','pos',0)
#particles.printPositions(LP_OutputDir+'/seededParticles','dat',0)
......
#Lagrangian particle-tracker.
#Requires simulationInfo.txt file with hydrodynamics simulation parameters:
# -> line1: simulation dt, line2: total nb of iters, line3: nb of iters between each export, line4: simulation start time, line5: output directory, line6: mesh name
#Mesh must be in planar coords in Meshes/ folder
#Makefile will download reef map into Data/ folder
#Hydrodynamics should be in Output/ folder
#Bathymetry (in msh format, same mesh as hydrodynamics) should be in Bath/ folder
from dgpy import *
import sys, time
import os
import os.path
from termcolor import colored
from ProjectMesh import *
print ''
print(colored('Lagrangian Advection-Diffusion particle tracking module for SLIM. Working offline.', "yellow"))
print ''
if (Msg.GetCommRank()==0):
print ''
print(colored('Compiling Libraries ....', "red"))
functionC.buildLibraryFromFile ("Libraries/GBR.cc", "Libraries/lib_gbr.so");
glib="Libraries/lib_gbr.so"
# Define simOutput directory & particle concentration
scratchFolder = 'Output'
simOutputDirName = ''
#simOutputDirName = 'Output200Kx3_4wks_windsOn_allTides_16112000'
simOutputDir = scratchFolder + simOutputDirName
LP_OutputDir = simOutputDir
concentration = 100.0 #No. of particles to seed per kmsq (100->1M particles)
reefMapFilename = "ReefMap_3600.dat"
mortOnOff = 1
XYZ = function.getCoordinates()
lonLatDegrees = functionC(glib,"lonLatDegrees",3,[XYZ])
print(colored('Loading simulation parameters ...',"red"))
simInfoFile = open(simOutputDir+'/simulationInfo.txt','r')
simInfo = simInfoFile.readlines()
simDt = float(simInfo[0])
simTotIter = int(simInfo[1])
simExportGap = int(simInfo[2])
simStartTime = int(simInfo[3])
#simOutputDir = simInfo[4].strip()
simMesh = 'gbr10K'#simInfo[5].strip()
simInfoFile.close()
#Lauch Makefile to generate meshes and download forcings
if(Msg.GetCommRank() == 0):
try : os.system("make "+reefMapFilename);
except: 0;
try : os.mkdir(LP_OutputDir);
except: 0;
simSolutions = [simExportGap*solutionNumber for solutionNumber in range(0,simTotIter/simExportGap+1)]
print(colored('Loading particle-tracking parameters ...',"red"))
LP_dt = 180. #(double) in seconds
LP_TotTime = 3*7*24.0 #(double) in hours
LP_Export = 10# int( (24.0*3600.0)/LP_dt ) #(int) export particle positions every ... steps
LP_Export_ConnectivityMatrix = int( (24.0*3600.0)/LP_dt ) #(int) export connectivity matrix every ... steps
LP_TotIter = 16# int(LP_TotTime*3600.0/LP_dt)
LP_StartTime = simStartTime # +24.0*3600#****+ (7*24.84*3600) #in seconds
LP_EndTime = LP_StartTime + LP_TotTime*3600.0
simItersPerLPIter = LP_dt/simDt
# Record simulation info to file
LP_SimInfo = ['LPT starts '+str((LP_StartTime-simStartTime)/3600)+' hrs after beginning of hydrodynamics simulation\n', 'LPT simulation length: '+str(LP_TotTime)+' hrs\n', 'Exporting every '+str(LP_Export)+' LPT steps\n', 'LPT dt: '+ str(LP_dt)+' s\n', 'No. of simulation iterations per LP iteration: '+str(simItersPerLPIter)+'\n']
LP_SimInfoFile = open(LP_OutputDir+'/LP_runInfo.txt','w')
LP_SimInfoFile.writelines(LP_SimInfo)
LP_SimInfoFile.close()
print(colored('Loading model ...',"red"))
model = GModel()
model.load("./Meshes/"+simMesh+"_tan.msh")
groups = dgGroupCollection(model, 2, 1) #model, dimension, order
#groups.buildGroupsOfInterfaces()
XYZ = function.getCoordinates()
bathDC = dgDofContainer(groups, 1);
if(os.path.exists("./Bath/"+simMesh+"_bath_smooth/"+simMesh+'_bath_smooth_COMP_0.msh')):
print(colored('Importing bathymetry.',"red"))
bathDC.importMsh("./Bath/"+simMesh+"_bath_smooth/"+simMesh+"_bath_smooth")
else:
print(colored('**** ERROR: Bathymetry file does not exist! ****',"red"))
print(colored('Defining diffusivity ...',"red"))
CCode = """
#include "fullMatrix.h"
#include "function.h"
#include "dgGroupOfElements.h"
#include "math.h"
#include "MElement.h"
extern "C" {
void diff(dataCacheMap *m, fullMatrix<double> &diff, fullMatrix<double> &alpha) {
double delta, K;
int ElementId;
const dgGroupOfElements *group = m->getGroupOfElements();
ElementId = m->getElementId();
MElement *e = group->getElement(ElementId);
delta = e->maxEdge();
K = alpha(0,0) * 0.00020551 * pow(delta,1.15);
for (int i=0; i<diff.size1(); i++) {
diff.set (i, 0, K);
}
}
}
"""
tmpCLib = "tempCLib.dylib"
if (Msg.GetCommRank() == 0 ) :
functionC.buildLibrary (CCode, tmpCLib)
alpha = functionConstant(1.0)
diffFC = functionC(tmpCLib,"diff",1,[alpha])
diffDC = dgDofContainer(groups,1)
diffDC.L2Projection(diffFC)
diffDC.exportMsh("./DIFFUSIVITY_RAW",0,0)
# Smooth diffusivity:
sys = linearSystemPETScDouble()
dofMan = dgDofManager.newCG (groups, 1, sys)
diffProjector = L2ProjectionContinuous(dofMan)
diffProjector.apply(diffDC,diffFC)
diffDC.exportMsh("./DIFFUSIVITY_smooooth",0,0)
# Load reef map.
print(colored('Loading reef map ...',"red"))
reefsDC = dgDofContainer(groups,1)
# 1. Create a reefs object (initialises object and creates reefs function)
reefsObj = reefs("./Data/ReefMap_Cropped_3600.dat", lonLatDegrees, 1)
# 2. Interpolate the function onto DC
reefsDC.interpolate(reefsObj)
reefsDC.exportMsh("./REEFS_RAW",0,0)
# 3. Delete reefs object data
reefsObj.clearReefsData()
dummyFM = fullMatrixDouble()
nbOfReefsFM = fullMatrixDouble()
reefsDC.minmax(dummyFM,nbOfReefsFM)
nbOfReefs = int(nbOfReefsFM(0,0)) + 1 #number of reefs = highest reef Id + 1 (lowest reefId = 0)
# Print size of smallest element:
maxEdge = 100000.0
maxMaxEdge = 0.0
minSurface = 1000000.0
maxSurface = 0.0
nbGroups = groups.getNbElementGroups()
for group in range( 0, nbGroups):
groupOfElements = groups.getElementGroup(group)
nbElements = groupOfElements.getNbElements()
for element in range( 0, nbElements):
el = groupOfElements.getElement(element)
if el.getVolume() < minSurface:
minSurface = el.getVolume()
if el.getVolume() > maxSurface:
maxSurface = el.getVolume()
if el.maxEdge() < maxEdge:
maxEdge = el.maxEdge()
if el.maxEdge() > maxMaxEdge:
maxMaxEdge = el.maxEdge()
print 'Smallest element\'s maxEdge is', maxEdge, 'm, Largest element\'s maxEdge is', maxMaxEdge, 'm'
print 'Smallest element\'s surface is', minSurface, 'msq, Largest element\'s surface is', maxSurface, 'msq'
print(colored('Seeding particles ...',"red"))
# ***SEED OVER ALL REEFS***
initCM = fullMatrixDouble( nbOfReefs, nbOfReefs )
particles = particleArray(groups, reefsDC, initCM, concentration)
particles.printPositions(LP_OutputDir+'/seededParticles','pos',0)
dummyPrinterICM = connectivityMatrixPrint()
dummyPrinterICM.CMPrinter(initCM, LP_OutputDir+'/connectivityMatrix_initial' )
#NB: clear dummy printer object
# ***SEED AT A POINT*** (NB: Look at method used to print results if changing seeding method)
#cartCoords = [0]*3
#geoCoordToCartesian(-19.9400,149.1100,cartCoords)
#projectPoint(cartCoords)
#particles = particleArray(groups, 1000, cartCoords[0], cartCoords[1], 0)
#particles.printPositions(LP_OutputDir+'/seededParticles','pos',0)
#particles.printPositions(LP_OutputDir+'/seededParticles','dat',0)
print ''
print '-----------------------'
print 'SIMULATION PARAMETERS:'
print 'Dt:', simDt,'s'
print 'Simulation length:', simTotIter, '(iter) |', simTotIter*simDt/3600.0, '(h)'
print 'Simulation starts:', simStartTime/3600.0, '(hr since origin) | Ends:', (simStartTime+simTotIter*simDt)/3600.0, '(hr since origin)'
print 'Data exported every', simExportGap, '(iter) |', simExportGap*simDt/60.0, '(mins)'
print '-----------------------'
print 'PARTICLE-TRACKING PARAMETERS:'
print 'Dt:', LP_dt, 's'
print 'Simulation length:', LP_TotIter, '(iter) |', LP_TotTime, '(h)'
print 'Simulation starts:', LP_StartTime/3600.0, '(hr since origin) | Ends:', LP_EndTime/3600.0, '(hr since origin)'
print 'Number of reefs in reef map:', nbOfReefs, '(highest reef Id on reefs DC)'
print 'Number of particles:', particles.printNbParticles()
print 'Mortality flag (1: ON, 0: OFF):', mortOnOff
LPsim_TimeOffset = LP_StartTime - simStartTime
print '-----------------------'
print 'Particle-tracker / Simulation time offset is:', LPsim_TimeOffset/3600.0
print '-----------------------'
print ''
simSolution = dgDofContainer(groups, 3)
simSolutionUnder = dgDofContainer(groups, 3)
simSolutionOver = dgDofContainer(groups, 3)
connectivityMatrix = fullMatrixDouble( nbOfReefs, nbOfReefs )
dummyPrinter = connectivityMatrixPrint()
#Create particle tracker object
print 'Initialising particle tracker ...'
particleTracker = dgParticleTracker2D(groups, particles, bathDC, diffDC, reefsDC, LP_TotIter, LP_dt, LP_OutputDir)
startcpu = time.clock()
for n in range(1,LP_TotIter+1):
t = simStartTime + LPsim_TimeOffset + float(n)*LP_dt
iterNumberWanted = int(LPsim_TimeOffset/simDt + n*simItersPerLPIter)
if iterNumberWanted in simSolutions:
simSolution.readMsh(simOutputDir+'/mr-gbr/mr-gbr-%06d'%( iterNumberWanted )+'.idx')
print '-> Found exact solution for iteration ', n,' (gbr-%06d'%( iterNumberWanted )+'). CPU time:', time.clock()-startcpu
else:
print '-> Calculating solution for n:', n, 'of',LP_TotIter+1,'| simIteration wanted:',iterNumberWanted, 'for t =',(t-simStartTime)/3600.0, '| CPU time:', time.clock()-startcpu
m=0
while (not iterNumberWanted+m in simSolutions):
m=m+1
# print 'Found simSolutionOver. m is:', m, '. File is (gbr-%06d'%( iterNumberWanted+m )+').'
t_solOver = simStartTime + float(iterNumberWanted+m)*simDt
simSolutionOver.setAll(10.0)
simSolutionOver.readMsh(simOutputDir+'/mr-gbr/mr-gbr-%06d'%( iterNumberWanted+m )+'.idx')
l=0
while (not iterNumberWanted-l in simSolutions):
l=l+1
# print 'Found simSolutionUnder. l is:', l, '. File is (gbr-%06d'%( iterNumberWanted-l )+').'
t_solUnder = simStartTime + float(iterNumberWanted-l)*simDt
simSolutionUnder.setAll(0.0)
simSolutionUnder.readMsh(simOutputDir+'/mr-gbr/mr-gbr-%06d'%( iterNumberWanted-l )+'.idx')
scaleFactor = 1.0 - (t-t_solUnder)/(t_solOver-t_solUnder)
# print 'Info: simSolutionUnder scale factor is ', scaleFactor
simSolutionUnder.scale(scaleFactor) #Scale simSolUnder (dof * relative dist. in t of sol. to simSolUnder)
simSolutionUnder.axpy(simSolutionOver,1.0-scaleFactor) #Scale simSolOver and add to scaled simSolUnder
simSolution.copy(simSolutionUnder) #Put solution in simSolution
#simSolution.exportMsh('./_LP-%06d'%(n), t, n)
#NB: 4th argument: Mortality (1: ON, 0: OFF)
particleTracker.particleMove(simSolution, connectivityMatrix, n, mortOnOff)
# Print output:
if ( n%(LP_Export) == 0 ):
particles.printPositions(LP_OutputDir+'/particleOutput_%06d'%( n ), 'pos', 0)
# *** More useful only for SEEDING AT SINGLE POINT ***
# if ( n%LP_Export == 0 ):
# particles.printPositions(LP_OutputDir+'/particleOutput_%06d'%( n ), 'dat', 0)
# particles.printPositionsFollowSingle(LP_OutputDir+'/particleOutputFollow', 'pos', n, LP_TotIter, LP_dt, 0)
# particles.printPositionsFollowSingle(LP_OutputDir+'/particleOutputFollow', 'pos', n, LP_TotIter, LP_dt, 400)
# particles.printPositionsFollowSingle(LP_OutputDir+'/particleOutputFollow', 'pos', n, LP_TotIter, LP_dt, 800)
if ( n%LP_Export_ConnectivityMatrix == 0 ):
dummyPrinter.CMPrinter(connectivityMatrix, LP_OutputDir+'/connectivityMatrix%06d'%( n ) )
print(colored('Particle tracking finished. Closing module.',"yellow"))
Msg.Exit(0)
/*
* General library of functions for GBR simulation.
* Part 1: "Numerical & computational" functions
* Part 2: "Physical" functions
*/
#include "fullMatrix.h"
#include "function.h"
#include "dgGroupOfElements.h"
#include "math.h"
#include "polynomialBasis.h"
static double pOx, pOy, pOz, pPhiX, pPhiY, pPhiZ, pThetaX, pThetaY, pThetaZ;
static bool init = false;
static void initialize() {
double phi = 146.5892/180*M_PI;
double theta = -17.147569/180*M_PI;
double R = 6.37101e6;
pOx = cos(theta)*cos(phi)*R;
pOy = cos(theta)*sin(phi)*R;
pOz = sin(theta)*R;
pPhiX = -sin(phi);
pPhiY = cos(phi);
pPhiZ = 0;
pThetaX = -sin(theta)*cos(phi);
pThetaY = -sin(theta)*sin(phi);
pThetaZ = cos(theta);
init = true;
}
extern "C" {
//Numerical & computational functions:
void lonLat (dataCacheMap *, fullMatrix<double> &lonlat, fullMatrix<double> &xyz) {
if (! init) initialize();
for (size_t i = 0; i< lonlat.size1(); i++) {
double x = pPhiX * xyz(i,0) + pThetaX * xyz(i,1) + pOx;
double y = pPhiY * xyz(i,0) + pThetaY * xyz(i,1) + pOy;
double z = pPhiZ * xyz(i,0) + pThetaZ * xyz(i,1) + pOz;
double r = sqrt(x*x+y*y+z*z);
lonlat.set(i, 0, atan2(y,x));
lonlat.set(i, 1, asin(z/r));
lonlat.set(i, 2, 0);
}
}
void latLon (dataCacheMap *, fullMatrix<double> &lonlat, fullMatrix<double> &xyz) {
if (! init) initialize();
for (size_t i = 0; i< lonlat.size1(); i++) {
double x = pPhiX * xyz(i,0) + pThetaX * xyz(i,1) + pOx;
double y = pPhiY * xyz(i,0) + pThetaY * xyz(i,1) + pOy;
double z = pPhiZ * xyz(i,0) + pThetaZ * xyz(i,1) + pOz;
double r = sqrt(x*x+y*y+z*z);
lonlat.set(i, 0, asin(z/r));
lonlat.set(i, 1, atan2(y,x));
lonlat.set(i, 2, 0);
}
}
void lonLatDegrees (dataCacheMap *, fullMatrix<double> &lonlat, fullMatrix<double> &xyz) {
if (! init) initialize();
for (size_t i = 0; i< lonlat.size1(); i++) {
double x = pPhiX * xyz(i,0) + pThetaX * xyz(i,1) + pOx;
double y = pPhiY * xyz(i,0) + pThetaY * xyz(i,1) + pOy;
double z = pPhiZ * xyz(i,0) + pThetaZ * xyz(i,1) + pOz;
double r = sqrt(x*x+y*y+z*z);
lonlat.set(i, 0, atan2(y,x)*180/M_PI);
lonlat.set(i, 1, asin(z/r)*180/M_PI);
lonlat.set(i, 2, 0);
}
}
void latLonDegrees (dataCacheMap *, fullMatrix<double> &lonlat, fullMatrix<double> &xyz) {
if (! init) initialize();
for (size_t i = 0; i< lonlat.size1(); i++) {
double x = pPhiX * xyz(i,0) + pThetaX * xyz(i,1) + pOx;
double y = pPhiY * xyz(i,0) + pThetaY * xyz(i,1) + pOy;
double z = pPhiZ * xyz(i,0) + pThetaZ * xyz(i,1) + pOz;
double r = sqrt(x*x+y*y+z*z);
lonlat.set(i, 1, atan2(y,x)*180/M_PI);
lonlat.set(i, 0, asin(z/r)*180/M_PI);
lonlat.set(i, 2, 0);
}
}
void lonLatVector(dataCacheMap *, fullMatrix<double> &uv, fullMatrix<double> &lonlat, fullMatrix<double> &u, fullMatrix<double> &v){
for (size_t i = 0; i< lonlat.size1(); i++) {
double lon=lonlat(i,0);
double lat=lonlat(i,1);
double uLon=u(i,0);
double uLat=v(i,0);
double lon0=146.5892*M_PI/180;
double lat0=-17.147569*M_PI/180;
double u1= -uLat*sin(lat)*cos(lon)-uLon*sin(lon);
double u2= -uLat*sin(lat)*sin(lon)+uLon*cos(lon);
double u3= uLat*cos(lat);
double alpha1= -sin(lon0);
double alpha2= cos(lon0);
double alpha3=0;
double beta1=-sin(lat0)*cos(lon0);
double beta2=-sin(lat0)*sin(lon0);
double beta3=cos(lat0);
uv.set(i, 0, u1*alpha1+u2*alpha2+u3*alpha3);
uv.set(i, 1, u1*beta1+u2*beta2+u3*beta3);
uv.set(i, 2, 0);
}
}
void latLonVector(dataCacheMap *, fullMatrix<double> &uv, fullMatrix<double> &latlon, fullMatrix<double> &u, fullMatrix<double> &v){
for (size_t i = 0; i< latlon.size1(); i++) {
double lat=latlon(i,0);
double lon=latlon(i,1);
double uLat=u(i,0);
double uLon=v(i,0);
double lon0=146.5892*M_PI/180;
double lat0=-17.147569*M_PI/180;
double u1= -uLat*sin(lat)*cos(lon)-uLon*sin(lon);
double u2= -uLat*sin(lat)*sin(lon)+uLon*cos(lon);
double u3= uLat*cos(lat);
double alpha1= -sin(lon0);
double alpha2= cos(lon0);
double alpha3=0;
double beta1=-sin(lat0)*cos(lon0);
double beta2=-sin(lat0)*sin(lon0);
double beta3=cos(lat0);
uv.set(i, 0, u1*alpha1+u2*alpha2+u3*alpha3);
uv.set(i, 1, u1*beta1+u2*beta2+u3*beta3);
uv.set(i, 2, 0);
}
}
void transport2velocity(dataCacheMap *,fullMatrix<double> &uv, fullMatrix<double> &UV, fullMatrix<double> &bath, fullMatrix<double> &solution){
for (size_t i = 0; i< solution.size1(); i++) {
uv.set(i,0,UV(i,0)/(bath(i,0)+solution(i,0)));
uv.set(i,1,UV(i,1)/(bath(i,0)+solution(i,0)));
uv.set(i,2,UV(i,2)/(bath(i,0)+solution(i,0)));
}
}
void computeEl(dataCacheMap *,fullMatrix<double> &val,fullMatrix<double> &solution) {
for (size_t i = 0; i< val.size1(); i++) {
val.set(i, 0, solution(i,0));
}
}
void computeVel(dataCacheMap *, fullMatrix<double> &val, fullMatrix<double> &solution){
for (size_t i = 0; i< val.size1(); i++) {
val.set(i, 0, solution(i,1));
val.set(i, 1, solution(i,2));
}
}
void computeErrorEl(dataCacheMap *,fullMatrix<double> &err,fullMatrix<double> &solution1, fullMatrix<double> &solution2) {
for (size_t i = 0; i< err.size1(); i++) {
err.set(i, 0, solution1(i,0)-solution2(i, 0));
}
}
void computeErrorVel(dataCacheMap *, fullMatrix<double> &err, fullMatrix<double> &solution1, fullMatrix<double> &solution2){
for (size_t i = 0; i< err.size1(); i++) {
err.set(i, 0, solution1(i,1)-solution2(i, 1));
err.set(i, 1, solution1(i,2)-solution2(i, 2));
}
}
void merge(dataCacheMap *,fullMatrix<double> &sol,fullMatrix<double> &eta, fullMatrix<double> &uv){
for (size_t i = 0; i< sol.size1(); i++) {
sol.set(i,0,eta(i,0));
sol.set(i,1,uv(i,0));
sol.set(i,2,uv(i,1));
}
}
void merge3(dataCacheMap *,fullMatrix<double> &sol,fullMatrix<double> &eta, fullMatrix<double> &uv, fullMatrix<double> &currentUV){
for (size_t i = 0; i< sol.size1(); i++) {
sol.set(i,0,eta(i,0));
sol.set(i,1,uv(i,0) + currentUV(i,1));
sol.set(i,2,uv(i,1) + currentUV(i,2));
}
}
//Physical functions:
void wind (dataCacheMap *, fullMatrix<double> &sol, fullMatrix<double> &xyz) {
for (size_t i = 0; i< sol.size1(); i++) {
sol.set(i,0,sin(xyz(i,1)/1e6)/1e6);
sol.set(i,1,0);
}
}
void coriolis (dataCacheMap *, fullMatrix<double> &sol, fullMatrix<double> &xyz) {
double theta = -17.147569/180*M_PI;
for (size_t i = 0; i< sol.size1(); i++) {
//sol.set(i,0,2*7.292e-5*xyz(i,2)/6.37101e6);
sol.set(i,0,2*7.292e-5*sin(theta));
}
}
void bottomDrag_old(dataCacheMap *, fullMatrix<double> &val, fullMatrix<double> &sol, fullMatrix<double> &bath){
for(int i=0; i<val.size1(); i++){
double H = sol(i, 0)+bath(i, 0);
double reef = 1.0;
if (H < 10.0) reef = 10.0;
val.set(i, 0, 9.81*reef*0.0235*0.0235/(pow(H, 1.333333333333)));
}
}
void bottomDrag(dataCacheMap *, fullMatrix<double> &val, fullMatrix<double> &sol, fullMatrix<double> &bath, fullMatrix<double> &reefs){
for(int i=0; i<val.size1(); i++){
double reef = 1.0;
double H = sol(i, 0)+bath(i, 0);
//Verify if there is a reef or not.
// Reef --> multiply normal bottom drag by factor 10x
if (reefs(i, 0) > 0.1) reef = 10.0;
val.set(i, 0, 9.81*reef*0.0235*0.0235/(pow(H, 1.333333333333)));
}
}
void smagorinsky(dataCacheMap *map, fullMatrix<double> &val, fullMatrix<double> &solGradient){
for(int i=0; i<val.size1(); i++){
double dudx = solGradient(i, 3);
double dudy = solGradient(i, 4);
double dvdx = solGradient(i, 6);
double dvdy = solGradient(i, 7);
const dgGroupOfElements *g = map->getGroupOfElements();
int elNum = map->getElementId();
double radi =g->getInnerRadius(elNum);
val(i, 0)= pow_int(0.1*radi, 2)*sqrt(2*dudx*dudx+2*dvdy*dvdy+pow_int(dudy+dvdx, 2));
}
}
void current (dataCacheMap *, fullMatrix<double> &current, fullMatrix<double> &solution) {
for (size_t i = 0; i< current.size1(); i++) {
current.set(i, 0, solution(i,1));
current.set(i, 1, solution(i,2));
current.set(i, 2, 0);
}
}
void windStress (dataCacheMap *,fullMatrix<double> &windstress,fullMatrix<double> &winduv, fullMatrix<double> &sol, fullMatrix<double> &bath) {
for (size_t i = 0; i< windstress.size1(); i++) {
double normWind=sqrt(winduv(i, 0)*winduv(i, 0)+winduv(i, 1)*winduv(i, 1));
double rho=1025;
double sol1=0.001*(0.630*normWind+0.066*normWind*normWind)*winduv(i, 0)/(rho*(sol(i, 0)+bath(i, 0)));
double sol2=0.001*(0.630*normWind+0.066*normWind*normWind)*winduv(i, 1)/(rho*(sol(i, 0)+bath(i, 0)));
windstress.set(i, 0, sol1);
windstress.set(i, 1, sol2);
}
}
void bathVec(dataCacheMap *, fullMatrix<double> &sol, fullMatrix<double> &bathymetry, fullMatrix<double> &normals, fullMatrix<double> &boundaryData) {
double vx, vy;
vx = boundaryData(0,1);
vy = boundaryData(0,2);
for (size_t i=0; i<bathymetry.size1(); i++) {
sol.set(i, 0, bathymetry(i,0) * (vx*normals(i,0) + vy*normals(i,1)) );
}
}
void SEC_UV(dataCacheMap *, fullMatrix<double> &UVout, fullMatrix<double> &bathInterfaceIntegrals, fullMatrix<double> &normals, fullMatrix<double> &boundaryData) {
double F, c, vx, vy;
F = boundaryData(0,0);
vx = boundaryData(0,1);
vy = boundaryData(0,2);
c = F / bathInterfaceIntegrals(0,0);
for (size_t i=0; i<UVout.size1(); i++) {
UVout.set(i, 0, 0); //Eta
UVout.set(i, 1, c * vx);
UVout.set(i, 2, c * vy);
}
}