Ulambators BoundaryBlock class (bndblock.cpp) v0.3
This class contains functions to create and store information of boundaries and interfaces
Contents:
-
Initialisation
-
Wall and Droplet creation
-
Droplet evolution
-
Special functions for drop merging, repulsion etc
-
Data writing
-
Drolet statistics
-
Remeshing procedures
-
Ubiquitous subroutines
The header optionally contains OpenMP support.
#include "bndblock.h" #include "bndelement.h" #ifdef _WIN32 #define _USE_MATH_DEFINES // for C++ #include <cmath> #endif #include <math.h> #include <iostream> #include <stdio.h> #include "combase.h" #include <stdlib.h> #include <fstream> #include <assert.h> #ifdef _WIN32 #include <Winsock2.h> #include <algorithm> #include <time.h> #endif using namespace std;Initialization
BndBlock::BndBlock(int blocks){ Elements = new FlowFace[blocks](); nbBlocks = blocks; runstep = 0; runtime = 0.0; nbPanels = -1; rlim = 0; //elementsPerLength; workdir = new char[120]; panel2block = new int[64]; // only 64 panels can exist, increase if mor are needed writeEnabled = 0; remesh_threshold = 2.0; SeedOption = 0; // 1 - Variable, else Constant IGP = new double[6]; IGW = new double[6]; IGP[0] = -0.932469514203152; IGP[1] = -0.661209386466265; IGP[2] = -0.238619186083197; IGP[3] = 0.238619186083197; IGP[4] = 0.661209386466265; IGP[5] = 0.932469514203152; IGW[0] = 0.171324492379170; IGW[1] = 0.360761573048139; IGW[2] = 0.467913934572691; IGW[3] = 0.467913934572691; IGW[4] = 0.360761573048139; IGW[5] = 0.171324492379170; } void BndBlock::Init(){ int counter = 0; ifid = 0; ifblock = 0; dynblock = 0; dynpanels = 0; nbPanels = 0; for (int k=0; k<nbBlocks; k++) { for (int j=nbPanels; j<(nbPanels+Elements[k].nbPanels); j++) { panel2block[j] = k; } nbPanels += Elements[k].nbPanels; } allBCX = (double*) calloc (nbPanels,sizeof(double)); allBCY = (double*) calloc (nbPanels,sizeof(double)); allBCZ = (double*) calloc (nbPanels,sizeof(double)); allXId = (double**) calloc (nbPanels,sizeof(double*)); allYId = (double**) calloc (nbPanels,sizeof(double*)); allBCtype = (int*) calloc (nbPanels,sizeof(int)); allSize = (int*) calloc (nbPanels,sizeof(int)); allRdist = (double**) calloc (nbPanels,sizeof(double*)); allCdist = (double**) calloc (nbPanels,sizeof(double*)); allWdist = (double**) calloc (nbPanels,sizeof(double*)); for (int k=0; k<nbBlocks; k++) { switch (Elements[k].BCkind) { case 0: // Outer Boundary for (int j=0; j<Elements[k].nbPanels; j++) { allBCX[counter] = Elements[k].BCvalueX[j]; allBCY[counter] = Elements[k].BCvalueY[j]; allBCtype[counter] = Elements[k].BCtype[j]; allXId[counter] = Elements[k].XId[j]; allYId[counter] = Elements[k].YId[j]; allSize[counter] = Elements[k].PanelSize[j]; counter++; } break; /* case 1: // Fibers for (int j=0; j<Elements[k].nbPanels; j++) { allBCX[counter] = 0.0; allBCY[counter] = 0.0; allBCZ[counter] = 0.0; allBCtype[counter] = 4; allXId[counter] = Elements[k].XId[j]; allYId[counter] = Elements[k].YId[j]; allSize[counter] = Elements[k].PanelSize[j]; counter++; } break; */ case 2: // Droplet allRdist[counter] = Elements[k].Rdist; allCdist[counter] = Elements[k].Cdist; allWdist[counter] = Elements[k].Wdist; for (int j=0; j<Elements[k].nbPanels; j++) { allBCX[counter] = 0.0; allBCY[counter] = 0.0; allBCtype[counter] = 21; if (Elements[k].BCtype[j]==22) { allBCtype[counter] = 22; } allXId[counter] = Elements[k].XId[j]; allYId[counter] = Elements[k].YId[j]; allSize[counter] = Elements[k].PanelSize[j]; counter++; } break; /* case 3: // Bordtrac allRdist[counter] = Elements[k].Rdist; allCdist[counter] = Elements[k].Cdist; allWdist[counter] = Elements[k].Wdist; for (int j=0; j<Elements[k].nbPanels; j++) { allBCX[counter] = 0.0; allBCY[counter] = 0.0; allBCtype[counter] = 23; allXId[counter] = Elements[k].XId[j]; allYId[counter] = Elements[k].YId[j]; allSize[counter] = Elements[k].PanelSize[j]; counter++; } break; */ default: cout<<"No valid boundary selected !"<<flush; exit(0); break; } if (Elements[k].BCkind<2) { ifid += Elements[k].nbPanels; ifblock++; } if (Elements[k].BCkind<1) { dynpanels += Elements[k].nbPanels; dynblock++; } } WrapAll(); allStock(); } void BndBlock::PanelInit(int blockId, int Size, int ofKind, int Btype, double *BCvalue, double *PosVec){ // One could check here that the number of panels is not exceeded int index; int numPanels = Elements[blockId].nbPanels; if (numPanels>=1) { index = Elements[blockId].PanelSize[numPanels-1]; if(!((fabs(Elements[blockId].XId[numPanels-1][index]-PosVec[0])<1e-6)&&(fabs(Elements[blockId].YId[numPanels-1][index]-PosVec[1])<1e-6))){ std::cout<<"ATTENTION: Geometry NOT closed, check boundary output!!!\n"; } } index = 0; for (int k = 0; k<numPanels; k++) { index += Elements[blockId].PanelSize[k]; } Elements[blockId].XId[numPanels] = &Elements[blockId].XE[index+5]; Elements[blockId].YId[numPanels] = &Elements[blockId].YE[index+5]; Elements[blockId].PanelSize[numPanels] = Size; switch (ofKind){ case 0: // Straight line for walls or in- outflow DrawLine(Elements[blockId].XId[numPanels], Elements[blockId].YId[numPanels], PosVec[0], PosVec[1], PosVec[2], PosVec[3], Size); break; case 1: // Curved line for walls Pos4 is curvature positive for outwards DrawSegm(Elements[blockId].XId[numPanels], Elements[blockId].YId[numPanels], PosVec[0], PosVec[1], PosVec[2], PosVec[3], PosVec[4], Size); break; case 2: // Pancake droplet DoCircle(Size, PosVec[2], Elements[blockId].XId[numPanels], Elements[blockId].YId[numPanels], PosVec[0], PosVec[1], int(PosVec[4]), PosVec[3]); break; case 3: // Elliptic droplet Pos0 is the major halfaxis and its inverse minor ha DoEllipse(Size, PosVec[2], PosVec[3], Elements[blockId].XId[numPanels], Elements[blockId].YId[numPanels], PosVec[0], PosVec[1]); break; case 4: // droplet with recangular middle section and caps, Pos0 is length of middle section DoEgg(Size, PosVec[2], PosVec[3], Elements[blockId].XId[numPanels], Elements[blockId].YId[numPanels], PosVec[0], PosVec[1]); break; case 5: // Droplet interface for attached boundaries DrawLine(&Elements[blockId].XId[numPanels][-1], &Elements[blockId].YId[numPanels][-1], PosVec[0], PosVec[1], PosVec[2], PosVec[3], Size+1); break; case 6: // Hyperbolic walls DrawHyperbel(Elements[blockId].XId[numPanels], Elements[blockId].YId[numPanels], PosVec[0], PosVec[1], PosVec[2], PosVec[3], Size); case 7: // Sinus shape walls with amplitude and wavenumber DrawSine(Elements[blockId].XId[numPanels], Elements[blockId].YId[numPanels], PosVec[0], PosVec[1], PosVec[2], PosVec[3], Size,PosVec[4], int(PosVec[5])); break; } Elements[blockId].BCtype[numPanels] = Btype; Elements[blockId].BCvalueX[numPanels] = BCvalue[0]; Elements[blockId].BCvalueY[numPanels] = BCvalue[1]; numPanels++; Elements[blockId].nbPanels = numPanels; if((numPanels>1)&&(numPanels==Elements[blockId].nbPanRes)){ if(!((fabs(Elements[blockId].XId[numPanels-1][Size]-Elements[blockId].XId[0][0])<1e-6)&&(fabs(Elements[blockId].YId[numPanels-1][Size]-Elements[blockId].YId[0][0])<1e-6))){ std::cout<<"ATTENTION: Geometry NOT closed, check boundary output!!!\n"; } Elements[blockId].XId[numPanels-1][Size+1] = Elements[blockId].XId[0][1]; Elements[blockId].YId[numPanels-1][Size+1] = Elements[blockId].YId[0][1]; } std::cout<<" Block "<<blockId<<" Panel "<<numPanels<<" of "<<Elements[blockId].nbPanRes<<" initialized.\n"<<flush; }Dedicated to droplet creation
void BndBlock::DoCircle(int prec, double radius, double *Xcir, double *Ycir, double XPos, double YPos, int Nmode, double Namp){Creates droplet, optional with perturbation single/multiple/noise
double radstep = 2*M_PI/prec; double radus; srand(time(NULL)); radus = 0.0; for(int n=0; n<=prec; n++){ if (Nmode>0) { radus = Namp*(cos(Nmode*n*radstep)); } else { radus = Namp*rand()/(float(RAND_MAX)+1); } Xcir[n] = radius*(1.0+radus)*cos(n*radstep)+ XPos; Ycir[n] = fabs(radius)*(1.0+radus)*sin(n*radstep)+ YPos; } } void BndBlock::DoEgg(int prec, double longeur, double largeur, double *Xcir, double *Ycir, double X0, double Y0){Make droplet which is elongated by a straight section – Option 4
double radius = largeur; double signum = 1.0-2.0*(longeur<0); longeur= fabs(longeur); int straight = (prec*longeur/(2*M_PI*radius+2*longeur)); int curved = prec/2-straight; double radstep = M_PI/curved; int n=0; int keep; for(int j=0; j<curved; j++){ Xcir[n] = -radius*cos(n*radstep)+X0; Ycir[n] = -signum*(longeur*1.0/2.0+radius*sin(n*radstep))+Y0; n++; } keep=n; for(int j=0; j<straight; j++){ Xcir[n] = radius+X0; Ycir[n] = signum*longeur*((n-keep)*1.0/straight-0.5)+Y0; n++; } keep = n; for(int j= 0; j<curved; j++){ Xcir[n] = radius*cos((n-keep)*radstep)+X0; Ycir[n] = signum*(longeur*1.0/2.0+radius*sin((n-keep)*radstep))+Y0; n++; } keep = n; for(int j=keep; j<prec; j++){ Xcir[n] = -radius+X0; Ycir[n] = signum*longeur*(0.5-(n-keep)*1.0/straight)+Y0; n++; } } void BndBlock::DoEllipse(int prec, double axisx, double axisy, double *Xcir, double *Ycir, double X0, double Y0){ // Make droplet which is elliptic major/minor axisx and axisy double radstep = 2*M_PI/prec; for(int n=0; n<prec; n++){ Xcir[n] = axisx*cos((n)*radstep)+X0; Ycir[n] = axisy*sin((n)*radstep)+Y0; } } void BndBlock::DrawLine(double *Xdraw, double *Ydraw, double X0, double Y0, double X1, double Y1, int nElm){ // Standard routine for straight lines refined at the end. Ovstp is tunable. double dx = (X1-X0); double dy = (Y1-Y0); double stlen = 2.0/nElm; double AQ, BQ, Ovstp,Sva; int keep=0; Ovstp = 1.1; if (Ovstp ==1) { double stlen = 1.0/nElm; for (int j=0; j<=nElm; j++) { Xdraw[j] = X0+dx*stlen*j; Ydraw[j] = Y0+dy*stlen*j; } } else { AQ = 2.5*(2-Ovstp)/(Ovstp-1); BQ = AQ+4/(2+AQ)-2; for (int j=0; j<nElm/2; j++) { Sva = j*1.0/nElm; // stlen = Sva; stlen = (AQ*Sva+1.0/(1.0+Sva*AQ)-1.0)/BQ; Xdraw[j] = X0+dx*stlen; Ydraw[j] = Y0+dy*stlen; keep = j; } for (int j=keep; j<=nElm; j++) { Sva = 1-j*1.0/nElm; // stlen = 1.0-Sva; stlen = 1.0 - (AQ*Sva+1.0/(1.0+Sva*AQ)-1.0)/BQ; Xdraw[j] = X0+dx*stlen; Ydraw[j] = Y0+dy*stlen; } } } void BndBlock::DrawSegm(double*Xdraw, double *Ydraw, double X0, double Y0, double X1, double Y1, double R, int nElm){ // Makes curved boundary segments, needs care when used double D = sqrt((X1-X0)*(X1-X0) + (Y1-Y0)*(Y1-Y0)); double Ym = -R/fabs(R)*sqrt(R*R - D*D/4); double gam = asin(D/(2*R)); double Xm = 0.5*(X0+X1) - Ym*(Y1-Y0)/D; double phi; double stlen = 1.0/nElm; double AQ, BQ, Ovstp,Sva; int keep=0; Ovstp = 1.3; AQ = 2.5*(2-Ovstp)/(Ovstp-1); BQ = AQ+4/(2+AQ)-2; Ym = 0.5*(Y0+Y1) + Ym*(X1-X0)/D; for (int j=0; j<nElm/2; j++) { Sva = j*1.0/nElm; stlen = (AQ*Sva+1.0/(1.0+Sva*AQ)-1.0)/BQ; phi = gam - stlen*gam*2.0; Xdraw[j] = Xm - R*(sin(phi)*(X1-X0)/D + cos(phi)*(Y1-Y0)/D); Ydraw[j] = Ym + R*(cos(phi)*(X1-X0)/D - sin(phi)*(Y1-Y0)/D); keep = j; } for (int j=keep; j<=nElm; j++) { Sva = 1-j*1.0/nElm; stlen = 1.0 - (AQ*Sva+1.0/(1.0+Sva*AQ)-1.0)/BQ; phi = gam - stlen*gam*2.0; Xdraw[j] = Xm - R*(sin(phi)*(X1-X0)/D + cos(phi)*(Y1-Y0)/D); Ydraw[j] = Ym + R*(cos(phi)*(X1-X0)/D - sin(phi)*(Y1-Y0)/D); } } void BndBlock::DrawSine(double*Xdraw, double *Ydraw, double X0, double Y0, double X1, double Y1, int nElm, double sinamp, int sinewave){ // Makes curved boundary segments, needs care when used double dx = (X1-X0); double dy = (Y1-Y0); double dabs = sqrt(dx*dx+dy*dy); double stlen = 2.0/nElm; double AQ, BQ, Ovstp,Sva; int keep=0; Ovstp = 1.0; if (Ovstp ==1) { double stlen = 1.0/nElm; for (int j=0; j<=nElm; j++) { Xdraw[j] = X0+dx*stlen*j+dy/dabs*sinamp*sin(sinewave*2.0*j*M_PI/nElm); Ydraw[j] = Y0+dy*stlen*j-dx/dabs*sinamp*sin(sinewave*2.0*j*M_PI/nElm); } } else { AQ = 2.5*(2-Ovstp)/(Ovstp-1); BQ = AQ+4/(2+AQ)-2; for (int j=0; j<nElm/2; j++) { Sva = j*1.0/nElm; // stlen = Sva; stlen = (AQ*Sva+1.0/(1.0+Sva*AQ)-1.0)/BQ; Xdraw[j] = X0+dx*stlen; Ydraw[j] = Y0+dy*stlen; keep = j; } for (int j=keep; j<=nElm; j++) { Sva = 1-j*1.0/nElm; // stlen = 1.0-Sva; stlen = 1.0 - (AQ*Sva+1.0/(1.0+Sva*AQ)-1.0)/BQ; Xdraw[j] = X0+dx*stlen; Ydraw[j] = Y0+dy*stlen; } } } void BndBlock::DrawHyperbel(double*Xdraw, double *Ydraw, double X0, double Y0, double X1, double Y1, int nElm){ double C = X0*Y0; assert(C==(Y1*X1)); double D1 = X0*X0-Y0*Y0; double D2 = X1*X1-Y1*Y1; double dd = (D2-D1)/nElm; for (int j=0; j<=nElm; j++){ Ydraw[j] = (( Y0 > 0 ) - ( Y0 < 0 ))*sqrt(sqrt(pow(D1+dd*j,2)/4+C*C)-(D1+dd*j)/2); Xdraw[j] = C/Ydraw[j]; } } void BndBlock::DropRead(int dropId, int timestep){ // Load droplet from solution used for visualization if(!writeEnabled){ std::cout<<"Writing is not enabled! No folder available.\nCall EnableWrite(dir) before ReadDrop(). Aborting program!\n"<<std::flush; exit(0); } char str[100]; char* thisdir = new char[120]; string kio, kio2; string delim = " "; // or " " double xpos, ypos; int npos, cutat; Elements[dropId].XId[0] = &Elements[dropId].XE[5]; Elements[dropId].YId[0] = &Elements[dropId].YE[5]; double *XId = Elements[dropId].XId[0]; double *YId = Elements[dropId].YId[0]; sprintf(thisdir,"%s/drop%its%i.dat", workdir.c_str(), dropId, timestep); fstream file(thisdir,ios::in); if (file.is_open()) { // Read Header file.getline(str,100); kio = str; cutat = kio.find_first_of(delim); kio2 = kio.substr(0,cutat); runtime = strtod(kio2.c_str(),NULL); kio = kio.substr(cutat+1); cutat = kio.find_first_of(delim); kio2 = kio.substr(0,cutat); Elements[dropId].Xrelativ = strtod(kio2.c_str(), NULL); kio = kio.substr(cutat+1); cutat = kio.find_first_of(delim); kio = kio.substr(0,cutat); Elements[0].Xrelativ = strtod(kio.c_str(), NULL); // Read Geometry npos =-1; while (file.eof()!=1) { npos++; file.getline(str,100); kio = str; cutat = kio.find_first_of(delim); kio2 = kio.substr(0,cutat); kio = kio.substr(cutat+1); cutat = kio.find_first_of(delim); kio = kio.substr(0,cutat); xpos = strtod(kio2.c_str(), NULL); ypos = strtod(kio.c_str(), NULL); XId[npos] = xpos; YId[npos] = ypos; } Elements[dropId].nbPanels = 1; Elements[dropId].PanelSize[0] = npos; runstep = timestep; if(nbPanels>=0){ allSize[ifid+dropId-ifblock] = Elements[dropId].PanelSize[0]; } file.close(); WrapInterface(dropId); Stock(dropId); std::cout<<"Drop Loaded with "<<npos<<" Elements.\n"; }else{ std::cout<<"Write directory could not be opened\n"; } }Dedicated to droplet evolution
void BndBlock::Evolve(int blockId, double Cold, double Cnew, double deltaT, double* xvec){ // Evolves droplets from buffer IntfXY to XYId int dynNodes; double *Xn, *Yn, *Xo, *Yo; dynNodes = 0; dynNodes = Elements[blockId].getFullSize(); Xn = Elements[blockId].XId[0]; Yn = Elements[blockId].YId[0]; Xo = &Elements[blockId].IntfX[5]; Yo = &Elements[blockId].IntfY[5]; for (int n=0; n<dynNodes; n++) { Xn[n] = Cold*Xo[n] + Cnew*Xn[n] + xvec[2*n]*deltaT;//*fabs(lnx); Yn[n] = Cold*Yo[n] + Cnew*Yn[n] + xvec[2*n+1]*deltaT;//*fabs(lny); } WrapInterface(blockId); } void BndBlock::Translate(int blockId, double Xdep, double Ydep, double Theta){ // Displace droplet by a solid body movement int dynNodes = Elements[blockId].getFullSize(); double* XA = Elements[blockId].XId[0]; double* YA = Elements[blockId].YId[0]; double* XT = &Elements[blockId].IntfX[5]; double* YT = &Elements[blockId].IntfY[5]; double BY, coz, xcm, ycm; getCoM(blockId, xcm, ycm); for (int n=0; n<dynNodes; n++) { BY = sqrt(pow(XT[n]-xcm,2)+pow(YT[n]-ycm,2)); coz = acos((XT[n]-xcm)/BY)*((YT[n]-ycm)>=0)+(2*M_PI-acos((XT[n]-xcm)/BY))*((YT[n]-ycm)<0); XA[n] = BY*cos(Theta+coz)+Xdep+xcm; YA[n] = BY*sin(Theta+coz)+Ydep+ycm; } WrapInterface(blockId); } void BndBlock::correctArea(int blockId){ // Corrects the Area by normal displacement of the interface double A1 = getArea(blockId); double A0 = Elements[blockId].initialArea; double circumference = 0; int dynNodes = Elements[blockId].getFullSize(); double* XA = Elements[blockId].XId[0]; double* YA = Elements[blockId].YId[0]; double* XT = &Elements[blockId].IntfX[5]; double* YT = &Elements[blockId].IntfY[5]; double lnx, lny, rr; if(A0!=0){ WrapInterface(blockId); Stock(blockId); for (int n=0; n<dynNodes; n++) { lnx = YA[n+1]-YA[n]; lny = XA[n+1]-XA[n]; rr = sqrt(lnx*lnx+lny*lny); circumference += rr; } circumference = (A0-A1)/circumference; for (int n=0; n<dynNodes; n++) { lnx = YT[n+1]-YT[n-1]; lny = XT[n-1]-XT[n+1]; rr = 1.0/sqrt(lnx*lnx+lny*lny); XA[n] += circumference*lnx*rr; YA[n] += circumference*lny*rr; } WrapInterface(blockId); } } void BndBlock::correctArea(){ for(int i=ifblock;i<nbBlocks;i++){ correctArea(i); } } void BndBlock::Stock(int blockId){ // Move droplet positins from XYId to buffer IntfXY int size = Elements[blockId].getFullSize(); for (int n=0; n<size+10;n++){ Elements[blockId].IntfX[n] = Elements[blockId].XE[n]; Elements[blockId].IntfY[n] = Elements[blockId].YE[n]; } } void BndBlock::UpdateCourbure(int blockId){ double *IX, *IY; int size; double dxl, dxr, dyl, dyr, delta1, delta2; IX = Elements[blockId].XId[0]; IY = Elements[blockId].YId[0]; size = Elements[blockId].getFullSize(); for(int n=0; n<size;n++){ dxl = IX[n]-IX[n-1]; dxr = IX[n+1]-IX[n]; dyl = IY[n]-IY[n-1]; dyr = IY[n+1]-IY[n]; delta1 = sqrt(dxl*dxl+dyl*dyl); delta2 = sqrt(dxr*dxr+dyr*dyr); dxl = -2.0*(dxr/delta2-dxl/delta1)/(delta1+delta2); dyl = -2.0*(dyr/delta2-dyl/delta1)/(delta1+delta2); Elements[blockId].KC[n] = sqrt(dxl*dxl+dyl*dyl)*(1.0-2.0*((dxl*dyr/delta2-dyl*dxr/delta2)<0)); } } void BndBlock::UpdateCourbure(){ for (int k=ifblock; k<nbBlocks; k++) { UpdateCourbure(k); } } void BndBlock::UpdateSize(){ int count=0; for (int n=ifblock; n<nbBlocks; n++) { for (int m=0; m<Elements[n].nbPanels; m++){ allSize[ifid+count] = Elements[n].PanelSize[m]; count++; } } } void BndBlock::WrapInterface(int blockId){ // Sets the ghost points int size = Elements[blockId].getFullSize(); double *Xp, *Yp; if(Elements[blockId].BCkind!=3){ for (int n=0; n<5; n++) { Elements[blockId].XId[0][n-5] = Elements[blockId].XId[0][size-5+n]; Elements[blockId].YId[0][n-5] = Elements[blockId].YId[0][size-5+n]; Elements[blockId].XId[0][size+n] = Elements[blockId].XId[0][n]; Elements[blockId].YId[0][size+n] = Elements[blockId].YId[0][n]; } } else { Xp = Elements[blockId].XId[0]; Yp = Elements[blockId].YId[0]; // Zero Angle - Problematic because implicit curvature assumes immobile ghostpoints Elements[blockId].XId[0][-2] = 2.0*Xp[-1]-Xp[0]; Elements[blockId].YId[0][-2] = 2.0*Yp[-1]-Yp[0]; Elements[blockId].YId[0][-3] = 3.0*Yp[-1]-2.0*Yp[0]; Elements[blockId].XId[0][-3] = 3.0*Xp[-1]-2.0*Xp[0]; Elements[blockId].XId[0][size+1] = 2.0*Xp[size]-Xp[size-1]; Elements[blockId].YId[0][size+1] = 2.0*Yp[size]-Yp[size-1]; Elements[blockId].YId[0][size+2] = 3.0*Yp[size]-2.0*Yp[size-1]; Elements[blockId].XId[0][size+2] = 3.0*Xp[size]-2.0*Xp[size-1]; } }Special functions for topography change currently disabled
void BndBlock::WallRepuls(){
double BndBlock::DropCorrect(int Drop1, int Drop2){
double BndBlock::DropApproach(int Drop1, int Drop2){
void BndBlock::PerformFusion(int Drop1, int Drop2, int proxL, int proxR){
int BndBlock::CheckBreakUp(int blockId, double kF){
void BndBlock::PerformCut(int blockId, int cutHere1, int cutHere2){
Dedicated to data writing
void BndBlock::EnableWrite(char* save2dir){ workdir = save2dir; writeEnabled = 1; } void BndBlock::WritePos(int blockId, int nbSlice){Write the droplet positions and curvature
if (writeEnabled) { char strL[128]; int nbEle; double* Xpos; double* Ypos; FILE* pfile; int saveshed = 0; if(blockId == 99){ saveshed = 1; blockId = 1; } nbEle = Elements[blockId].getFullSize(); Xpos = Elements[blockId].XId[0]; Ypos = Elements[blockId].YId[0]; // std::cout<<" Elements present "<<nbEle<<" "; if (saveshed) { sprintf(strL,"%s/drop1ts999.dat",workdir.c_str()); pfile = fopen(strL,"w"); UpdateCourbure(blockId); fprintf(pfile, "%1.15f %1.15f %1.15f\n", runtime, Elements[blockId].Xrelativ, Elements[blockId].Yrelativ); for (int n=nbSlice; n<nbEle; n++) { fprintf(pfile, "%1.15f %1.15f %1.15f\n",Xpos[n],Ypos[n],Elements[blockId].KC[n]); } }else if (Elements[blockId].BCkind<1) { sprintf(strL,"%s/bnd%ipos.dat", workdir.c_str(), blockId); pfile = fopen(strL,"w"); for (int j=0; j<nbEle; j++) { fprintf(pfile, "%1.15f %1.15f \n",Xpos[j],Ypos[j]); } } else { sprintf(strL,"%s/drop%its%d.dat",workdir.c_str(), blockId, nbSlice); pfile = fopen(strL,"w"); UpdateCourbure(blockId); fprintf(pfile, "%1.15f %1.15f %1.15f\n", runtime, Elements[blockId].Xrelativ, Elements[blockId].Yrelativ); for (int n=0; n<nbEle; n++) { fprintf(pfile, "%1.15f %1.15f %1.15f\n",Xpos[n],Ypos[n],Elements[blockId].KC[n]); } } fclose(pfile); cout <<"Position "<<nbSlice<<" written for "<<blockId<<" object"<<" at time "<<runtime<<"\n"; } } void BndBlock::LogAreaPos(int logId, int blockId){ char logurl[120]; FILE* logfile; double A0, Xp, Yp; A0 = getCoM(blockId, Xp, Yp); sprintf(logurl,"%s/log%i.txt", workdir.c_str(),logId); logfile = fopen(logurl, "a"); fprintf(logfile, "%2.9f %2.9f %2.9f\n", A0, Xp, Yp); fclose(logfile); } void BndBlock::LogTimeStep(int logId){ char logurl[120]; FILE* logfile; sprintf(logurl,"%s/log%i.txt", workdir.c_str(),logId); logfile = fopen(logurl, "a"); fprintf(logfile, "%i %f ", runstep, runtime); fclose(logfile); } double BndBlock::LogCustom(int logid, int modusid){To write a log file with customized data. Delete or add routines as you like.
char logurl[120]; FILE* logfile; double xcm, ycm, perimeter; int imax; switch(modusid){ case 0: break; case 1: { // FOR ARTULAM double perimeter = 0; for(int i=0; i<Elements[1].PanelSize[0];i++){ perimeter += sqrt( pow(Elements[1].XId[0][i+1]-Elements[1].XId[0][i] ,2)+pow(Elements[1].YId[0][i+1]-Elements[1].YId[0][i],2)); } sprintf(logurl,"%s/log%i.txt", workdir.c_str(),logid); logfile = fopen(logurl, "a"); getCoM(1, xcm, ycm); someWay += xcm; fprintf(logfile, "%f %f %f %f\n", xcm, ycm, perimeter, runtime); fclose(logfile); return xcm; break; } case 2: { // FOR CORDERO STRETCH perimeter = 0; double xmax = -10.0; double xmin = 10.0; double ymax = -10.0; double ymin = 10.0; for(int i=0; i<Elements[1].PanelSize[0];i++){ perimeter += sqrt( pow(Elements[1].XId[0][i+1]-Elements[1].XId[0][i] ,2)+pow(Elements[1].YId[0][i+1]-Elements[1].YId[0][i],2)); xmax = max(xmax, Elements[1].XId[0][i]); xmin = min(xmin, Elements[1].XId[0][i]); ymax = max(ymax, Elements[1].YId[0][i]); ymin = min(ymin, Elements[1].YId[0][i]); } char logurl[120]; FILE* logfile; sprintf(logurl,"%s/log%i.txt", workdir.c_str(),logid); logfile = fopen(logurl, "a"); getCoM(1, xcm, ycm); someWay += xcm; fprintf(logfile, "%f %f %f %f %f %f\n", xcm, ycm, xmax-xmin, ymax-ymin, perimeter, runtime); fclose(logfile); return xcm; break; } case 3: { //get extreme points for relax study double ynord, ysud, xest, xouest, kapnord, kapsud, kapest, kapouest; int blockid = ifblock; // Droplet Element sprintf(logurl,"%s/log%i.txt", workdir.c_str(),logid); logfile = fopen(logurl, "a"); double nordmin, sudmin, estmin, ouestmin; int imin, jmax, jmin; double *XA = Elements[blockid].XId[0]; double *YA = Elements[blockid].YId[0]; nordmin = 10; sudmin = 10; estmin =10; ouestmin = 10; for (int j=0; j<Elements[blockid].PanelSize[0]; j++) { if ((YA[j]>0)&&(fabs(XA[j])<nordmin)) { nordmin = fabs(XA[j]); jmax = j; } if ((YA[j]<0)&&(fabs(XA[j])<sudmin)) { sudmin = fabs(XA[j]); jmin = j; } if ((XA[j]>0)&&(fabs(YA[j])<estmin)) { estmin = fabs(YA[j]); imax = j; } if ((XA[j]<0)&&(fabs(YA[j])<ouestmin)) { ouestmin = fabs(YA[j]); imin = j; } } UpdateCourbure(blockid); kapnord = Elements[blockid].KC[jmax]; kapsud = Elements[blockid].KC[jmin]; kapest = Elements[blockid].KC[imax]; kapouest = Elements[blockid].KC[imin]; ynord = YA[jmax]; ysud = YA[jmin]; xest = XA[imax]; xouest = XA[imin]; fprintf(logfile, "%f %f %f %f %f %f %f %f %f\n", runtime, kapnord, kapsud, kapest, kapouest, ynord, ysud, xest, xouest); fclose(logfile); return kapest; break; } } return 0; } void BndBlock::DisableWrite(){ writeEnabled = 0; }Dedicated to determine basic droplet statistics
double BndBlock::getArea(int blockId){Integration of the droplet Area
double xgs, ygs, dx, dy, r, lnx, lny; int nbElms = Elements[blockId].getFullSize(); double A = 0.0; double *Xc = Elements[blockId].XId[0]; double *Yc = Elements[blockId].YId[0]; if (Elements[blockId].BCkind!=3) { WrapInterface(blockId); for(int n=0; n<nbElms; n++){ xgs = 0.5*(Xc[n+1]+Xc[n]); ygs = 0.5*(Yc[n+1]+Yc[n]); dx = Xc[n+1]-Xc[n]; dy = Yc[n+1]-Yc[n]; r = dx*dx+dy*dy; r = sqrt(r); lnx = dy/r; lny = -dx/r; A += r*(xgs*lnx+ygs*lny); } A = A*0.5; // because div(x) = 2 and GW6 integrates over 2 length }else{ for(int n=-1; n<=nbElms; n++){ xgs = 0.5*(Xc[n+1]+Xc[n]); ygs = 0.5*(Yc[n+1]+Yc[n]); dx = Xc[n+1]-Xc[n]; dy = Yc[n+1]-Yc[n]; r = dx*dx+dy*dy; r = sqrt(r); lnx = dy/r; lny = -dx/r; A += (ygs*lny+xgs*lnx)*r; } A = A*0.5; } return A; } double BndBlock::getCoM(int blockId, double& Xcm, double& Ycm){Integration of center of area
double xgs, ygs, dx, dy, r, lnx, lny; int nbElms; double* Xc = Elements[blockId].XId[0]; double* Yc = Elements[blockId].YId[0]; double A = 0.0; nbElms = Elements[blockId].getFullSize(); A = 0; Xcm = 0; Ycm = 0; for(int n=0; n<nbElms; n++){ xgs = 0.5*(Xc[n+1]+Xc[n]); ygs = 0.5*(Yc[n+1]+Yc[n]); dx = Xc[n+1]-Xc[n]; dy = Yc[n+1]-Yc[n]; r = dx*dx+dy*dy; r = sqrt(r); lnx = dy/r; lny = -dx/r; A += ygs*lny*r; A -= xgs*lnx*r; Xcm += r*xgs*ygs*lny; Ycm += r*0.5*ygs*ygs*lny; } Xcm = Xcm/A; Ycm = Ycm/A; return A; } double BndBlock::getHiMo(int blockId, double& Xcm, double& Ycm, double& Mcm){Integration of center of area
double xgs, ygs, dx, dy, x1, y1, r, lnx, lny; int nbElms; double* Xc = Elements[blockId].XId[0]; double* Yc = Elements[blockId].YId[0]; double A = 0.0; Xcm = 0.0; Ycm = 0.0; Mcm = 0.0; nbElms = Elements[blockId].getFullSize(); for(int n=0; n<nbElms; n++){ xgs = 0.5*(Xc[n+1]+Xc[n]); ygs = 0.5*(Yc[n+1]+Yc[n]); dx = Xc[n+1]-Xc[n]; dy = Yc[n+1]-Yc[n]; r = dx*dx+dy*dy; r = sqrt(r); lnx = dy/r; lny = -dx/r; for (int i=0; i<6; i++) { x1 = xgs+dx*IGP[i]*0.5; y1 = ygs+dy*IGP[i]*0.5; A += r*(x1*lnx+y1*lny)*0.5*IGW[i]; Xcm += r*(y1*x1*lny)*0.5*IGW[i]; Ycm += r*(y1*x1*lnx)*0.5*IGW[i]; // Mcm += r*(x1*y1*y1*lnx+x1*x1*y1*lny); } } Xcm = 2.0*Xcm/A; Ycm = 2.0*Ycm/A; A = -A/2.0; for(int n=0; n<nbElms; n++){ xgs = 0.5*(Xc[n+1]+Xc[n])-Xcm; ygs = 0.5*(Yc[n+1]+Yc[n])-Ycm; dx = Xc[n+1]-Xc[n]; dy = Yc[n+1]-Yc[n]; r = dx*dx+dy*dy; r = sqrt(r); lnx = dy/r; lny = -dx/r; for (int i=0; i<6; i++) { x1 = xgs+dx*IGP[i]*0.5; y1 = ygs+dy*IGP[i]*0.5; Mcm += r*(x1*y1*y1*lnx+x1*x1*y1*lny)*IGW[i]*0.5; } } Mcm = -Mcm; return A; } void BndBlock::LoadNOrder(int blockId){Compute 5th order polynomial of the droplet interfaces No longer functional because no memory assigned to polynomial coefficients
int order = 5; double *FaceX; double *FaceY; int size; double *AR = new double[order*order]; double *rhsx = new double[order]; double *rhsy = new double[order]; double dx,dy,rp; int ip,rs, frac; FaceX = Elements[blockId].XId[0]; FaceY = Elements[blockId].YId[0]; size = Elements[blockId].getFullSize(); for (int k = 0; k<size; k++) { frac = (order+0.5)/2; ip = 0; rs = -1.0; for (int i = 0; i<order; i++) { if(i==frac){ ip++; rs = 1; } dx = FaceX[k+ip-frac]-FaceX[k]; dy = FaceY[k+ip-frac]-FaceY[k]; rp = rs*sqrt(dx*dx+dy*dy); rhsx[i] = dx; rhsy[i] = dy; for (int j = 0; j<order; j++) { AR[i+j*order] = pow(rp,j+1); } ip++; } dgsev2rhs(AR, rhsx, rhsy, order); } }Dedicated to Remeshing
void BndBlock::SeedDrop(){ if (rlim>0.0) { for (int j=ifblock; j<nbBlocks; j++) { SeedEquiDist(j, rlim); } } } void BndBlock::RemeshOn(double targetsize, int modus){ rlim = targetsize; SeedOption = modus; } void BndBlock::SeedEquiDist(int blockId, double limp){Redistributes the droplet elements uniformly
int order = 3; double *XA, *XB; double *YA, *YB; double* AR = new double[order*order]; double *ace = new double[order]; double *bce = new double[order]; double rd, dx, dy, newx, newy, rmin, rnorm,rp; int jump, size, ReDist; int decal = 0; double rmax; double dsremain = 0; int nremain=0; double difadap=0; int ip,rs, frac; XB = Elements[blockId].XId[0]; YB = Elements[blockId].YId[0]; XA = &Elements[blockId].IntfX[5]; YA = &Elements[blockId].IntfY[5]; jump = 0; size = Elements[blockId].getFullSize(); ReDist = 0; rmax = remesh_threshold*limp; rmin = limp/remesh_threshold; dx = XB[-1]-XB[0]; dy = YB[-1]-YB[0]; rd = sqrt(dx*dx+dy*dy); if((rd>rmax)||(rd<rmin)){ ReDist = 1; } int kstart = 1; if(Elements[blockId].BCkind==3){ jump--; kstart=0; } for (int k=kstart; k<=size; k++) { // k should be 1 dx = XB[k]-XB[k-1]; dy = YB[k]-YB[k-1]; rd = sqrt(dx*dx+dy*dy); if((rd>rmax)||(rd<rmin)){ ReDist = 1; break; } } if(ReDist){ Stock(blockId); // LoadNOrder(blockId); //Copy Point XB[decal] = XA[kstart-1]; YB[decal] = YA[kstart-1]; rnorm = limp; for (int k=kstart; k<=(size); k++) { // k should be 1 dx = XA[k]-XA[k-1]; dy = YA[k]-YA[k-1]; rd = sqrt(dx*dx+dy*dy);This piece of code smoothes over the last elements the distance so no reconnection no size contrast exists
if (k>=size-10) { dsremain = limp-rnorm; for (int j=k; j<=size; j++) { dx = XA[k]-XA[k-1]; dy = YA[k]-YA[k-1]; dsremain += sqrt(dx*dx+dy*dy); } nremain = int(round(dsremain/limp)); difadap = (dsremain/(nremain*limp)); rnorm += (difadap-1.0)*limp; limp = limp*difadap; } { frac = int((order+0.5)/2.0); ip = 0; rs = -1.0; for (int i = 0; i<order; i++) { if(i==frac){ ip++; rs = 1; } dx = XA[k+ip-frac-1]-XA[k-1]; dy = YA[k+ip-frac-1]-YA[k-1]; rp = rs*sqrt(dx*dx+dy*dy); ace[i] = dx; bce[i] = dy; for (int j = 0; j<order; j++) { AR[i+j*order] = pow(rp,j+1); } ip++; } } dgsev2rhs(AR, ace, bce, order); while (rnorm<rd){ newx = XA[k-1]; newy = YA[k-1]; for (int j=0; j<order; j++) { newx += ace[j]*pow(rnorm, j+1); newy += bce[j]*pow(rnorm, j+1); } jump++; XB[jump+decal] = newx; YB[jump+decal] = newy; rnorm += limp; } rnorm += -rd; } if (Elements[blockId].BCkind==3) { rd = sqrt(pow(XB[jump]-XA[size],2) + pow(YB[jump]-YA[size],2)); } else { rd = sqrt(pow(XB[jump+decal]-XB[decal],2) + pow(YB[jump+decal]-YB[decal],2)); } if(Elements[blockId].BCkind!=3){ jump += (rd>(0.75*limp)); for (int k=0; k<decal; k++) { XB[k] = XB[jump+k]; YB[k] = YB[jump+k]; } } else { jump += (rd>(0.5*limp)); XB[-1] = XA[-1]; YB[-1] = YA[-1]; XB[jump] = XA[size]; YB[jump] = YA[size]; } Elements[blockId].PanelSize[0] = jump; std::cout<<"\br"; } allSize[ifid+blockId-ifblock] = Elements[blockId].PanelSize[0]; WrapInterface(blockId); Stock(blockId); }void BndBlock::SeedVarDist(int blockId, double rliming){
Ubiquitous subroutines
void BndBlock::allWrite(){ for (int k=0; k<nbBlocks; k++) { if((Elements[k].BCkind>0)||(runstep==0)){ WritePos(k, runstep); } } } void BndBlock::allStock(){ for (int k=0; k<nbBlocks; k++) { if(Elements[k].BCkind>0){ Stock(k); } } } void BndBlock::allEvolve(double weightOldPos, double weightNewPos, double deltaT, double *Uvec){ int index = 0; /* for (int k=0; k<ifblock; k++) { index += 2*Elements[k].getFullSize(); } */ for (int k=ifblock; k<nbBlocks; k++) { Evolve(k, weightOldPos, weightNewPos, deltaT, &Uvec[index]); index += 2*Elements[k].getFullSize(); } } void BndBlock::WrapAll(){ for (int k=dynblock; k<nbBlocks; k++) { WrapInterface(k); } } int BndBlock::getFixSize(){ int size=0; for (int k=0; k<nbBlocks; k++) { if(Elements[k].BCkind>0){ break; } size += Elements[k].getFullSize(); } return size; } int BndBlock::getAllSize(){ int size=0; for (int k=0; k<nbBlocks; k++) { size += 2*Elements[k].getFullSize(); if (Elements[k].BCkind==1) { size += 3; } } return size; }Created by Mathias Nagel on 31.05.12.
Copyright (c) 2012 EPFL. All rights reserved.