Skip to content

GitLab

Commit b000b0d4 authored by Phillip Lino Rall's avatar Phillip Lino Rall
Browse files

finished merging christops with phillips version : grid_transform_phillip is up to date now

parent 31464f21
Hide whitespace changes
Inline Side-by-side
program/source/grid/blockgrid.cc
View file @ b000b0d4
......@@ -42,25 +42,6 @@
int Blockgrid::id_count_grid = 0;
bool compareVectors(D3vector v, D3vector w)
{
if (!(v == w))
{
cout << "should: ";
v.Print();
cout<<endl;
cout << "but is: ";
w.Print();
cout<<endl;
cout<<endl;
return true;
}
else
{
// cout << "fine " << endl;
return false;
}
}
Blockgrid::Blockgrid() {
id_of_grid = id_count_grid; ++id_count_grid;
......@@ -116,15 +97,6 @@ Blockgrid::Blockgrid ( Unstructured_grid *ug_, int Nx, int Ny, int Nz ) {
number_points[i] = Nz;
}
// Give_Nx_hexahedron(0);
// blockgrid_Coordinates.resize(ug->Give_number_hexahedra());
// blockgrid_Coordinates_set.resize(ug->Give_number_hexahedra());
// cout << "still test here!!!" << endl;
// for (int iter = 0 ; iter <ug->Give_number_hexahedra() ; iter++ )
// {
// blockgrid_Coordinates.at(iter).resize(Give_Nx_hexahedron(iter)*Give_Ny_hexahedron(iter)*Give_Nz_hexahedron(iter));
// blockgrid_Coordinates_set.at(iter).resize(Give_Nx_hexahedron(iter)*Give_Ny_hexahedron(iter)*Give_Nz_hexahedron(iter));
// }
variable_set = false;
}
......@@ -543,7 +515,6 @@ D3vector Blockgrid::Give_coord_hexahedron ( int id_hex,
D3vector QPsi[6];
//if (transformFromQuadrangle)
if(ug->Give_transform_From_Quadrangle())
{
//geht bisher nur fuer 1 hexaeder, allgemien so richtig?
......
program/source/grid/blockgrid.h
View file @ b000b0d4
......@@ -147,7 +147,6 @@ class Blockgrid {
bool transformFromQuadrangle;
bool variable_set; // set true falls eine Variable konstruiert
int id_of_grid;
......
program/source/interpol/interpol.cc
View file @ b000b0d4
......@@ -30,35 +30,14 @@
#include "../extemp/extemp.h"
#include "../extemp/parallel.h"
#include "../extemp/variable.h"
#include "../mympi.h"
#include "../abbrevi.h"
#include "../parameter.h"
#include "../math_lib/math_lib.h"
#include "../math_lib/mg_array.h"
#include "../basics/basic.h"
#include "../grid/elements.h"
#include "../grid/parti.h"
#include "../grid/ug.h"
#include "../grid/examples_ug.h"
#include "../grid/blockgrid.h"
#include "../grid/marker.h"
#include "../extemp/extemp.h"
#include "../extemp/parallel.h"
#include "../extemp/variable.h"
#include "../extemp/cellvar.h"
#include "../extemp/co_fu.h"
#include "../extemp/functor.h"
#include "interpol.h"
#include <sstream>
#include <iomanip>
#include "assert.h"
/////////////////////////////////////////////////////////////
// 1. Interpolate from blockgrid to rectangular blockgrid
/////////////////////////////////////////////////////////////
......@@ -121,24 +100,33 @@ Intermadiate_grid_for_PointInterpolator::Intermadiate_grid_for_PointInterpolator
}
*/
Interpolate_on_structured_grid::~Interpolate_on_structured_grid() {
delete[] ids_hex;
delete[] ids_i;
delete[] ids_j;
delete[] ids_k;
delete[] typ_tet;
delete[] lambda;
}
Interpolate_on_structured_grid::Interpolate_on_structured_grid(int nx_, int ny_, int nz_,
D3vector pWSD, D3vector pENT,
Blockgrid& blockgrid_) {
D3vector pWSD, D3vector pENT,
Blockgrid& blockgrid_) {
int Nx, Ny, Nz;
int typ;
// int typ;
assert(nx_ > 1);
assert(ny_ > 1);
assert(nz_ > 1);
int ilmin, jlmin, klmin;
int ilmax, jlmax, klmax;
//int ilmin, jlmin, klmin;
//int ilmax, jlmax, klmax;
double factor = 0.1;
// double factor = 0.00001;
D3vector lam;
// D3vector lam;
blockgrid = &blockgrid_;
ug = blockgrid->Give_unstructured_grid();
......@@ -149,28 +137,30 @@ Interpolate_on_structured_grid::Interpolate_on_structured_grid(int nx_, int ny_,
if(nx_>1)
hx = (pENT.x - pWSD.x) / (nx_-1);
else
else
hx = 1.0;
if(ny_>1)
hy = (pENT.y - pWSD.y) / (ny_-1);
else
else
hy = 1.0;
if(nz_>1)
hz = (pENT.z - pWSD.z) / (nz_-1);
else
else
hz = 1.0;
int num_total = nx * ny * nz;
/*
D3vector cWSD, cESD;
D3vector cWND, cEND;
D3vector cWST, cEST;
D3vector cWNT, cENT;
*/
D3vector boxWSD, boxENT;
// D3vector boxWSD, boxENT;
D3vector ploc;
// D3vector ploc;
ids_hex = new int[num_total];
......@@ -189,59 +179,61 @@ Interpolate_on_structured_grid::Interpolate_on_structured_grid(int nx_, int ny_,
Ny = blockgrid->Give_Ny_hexahedron(id_hex);
Nz = blockgrid->Give_Nz_hexahedron(id_hex);
#pragma omp parallel for num_threads(UGBlocks::numThreadsToTake) if(UGBlocks::useOpenMP)
for(int k=0;k<Nz;++k)
for(int j=0;j<Ny;++j)
for(int i=0;i<Nx;++i) {
for(int j=0;j<Ny;++j)
for(int i=0;i<Nx;++i) {
// corner points of general hex-cell
cWSD = blockgrid->Give_coord_hexahedron(id_hex,i, j, k );
cESD = blockgrid->Give_coord_hexahedron(id_hex,i+1,j ,k );
cWND = blockgrid->Give_coord_hexahedron(id_hex,i, j+1,k );
cEND = blockgrid->Give_coord_hexahedron(id_hex,i+1,j+1,k );
cWST = blockgrid->Give_coord_hexahedron(id_hex,i, j, k+1);
cEST = blockgrid->Give_coord_hexahedron(id_hex,i+1,j ,k+1);
cWNT = blockgrid->Give_coord_hexahedron(id_hex,i, j+1,k+1);
cENT = blockgrid->Give_coord_hexahedron(id_hex,i+1,j+1,k+1);
// bounding box calculation
boxWSD.x = MIN(MIN(MIN(cWSD.x,cESD.x),MIN(cWND.x,cEND.x)),
MIN(MIN(cWST.x,cEST.x),MIN(cWNT.x,cENT.x))) - factor *hx;
boxWSD.y = MIN(MIN(MIN(cWSD.y,cESD.y),MIN(cWND.y,cEND.y)),
MIN(MIN(cWST.y,cEST.y),MIN(cWNT.y,cENT.y))) - factor *hy;
boxWSD.z = MIN(MIN(MIN(cWSD.z,cESD.z),MIN(cWND.z,cEND.z)),
MIN(MIN(cWST.z,cEST.z),MIN(cWNT.z,cENT.z))) - factor *hz;
boxENT.x = MAX(MAX(MAX(cWSD.x,cESD.x),MAX(cWND.x,cEND.x)),
MAX(MAX(cWST.x,cEST.x),MAX(cWNT.x,cENT.x))) + factor *hx;
boxENT.y = MAX(MAX(MAX(cWSD.y,cESD.y),MAX(cWND.y,cEND.y)),
MAX(MAX(cWST.y,cEST.y),MAX(cWNT.y,cENT.y))) + factor *hy;
boxENT.z = MAX(MAX(MAX(cWSD.z,cESD.z),MAX(cWND.z,cEND.z)),
MAX(MAX(cWST.z,cEST.z),MAX(cWNT.z,cENT.z))) + factor *hz;
// calculation of indices of a collection of cells of structured grid which contains bounding box
ilmin = Ganzzahliger_Anteil((boxWSD.x - pWSD.x) / hx);
jlmin = Ganzzahliger_Anteil((boxWSD.y - pWSD.y) / hy);
klmin = Ganzzahliger_Anteil((boxWSD.z - pWSD.z) / hz);
ilmax = Ganzzahliger_Anteil((boxENT.x - pWSD.x) / hx);
jlmax = Ganzzahliger_Anteil((boxENT.y - pWSD.y) / hy);
klmax = Ganzzahliger_Anteil((boxENT.z - pWSD.z) / hz);
/*
cout << " indices: "
<< " ilmin: " << ilmin
<< " jlmin: " << jlmin
<< " klmin: " << klmin
<< " ilmax: " << ilmax
<< " jlmax: " << jlmax
<< " klmax: " << klmax
<< " boxWSD.x: " << boxWSD.x
<< " cWSD.x: " << cWSD.x
<< " Nx: " << Nx
<< endl;
*/
/*
D3vector cWSD = blockgrid->Give_coord_hexahedron(id_hex,i, j, k );
D3vector cESD = blockgrid->Give_coord_hexahedron(id_hex,i+1,j ,k );
D3vector cWND = blockgrid->Give_coord_hexahedron(id_hex,i, j+1,k );
D3vector cEND = blockgrid->Give_coord_hexahedron(id_hex,i+1,j+1,k );
D3vector cWST = blockgrid->Give_coord_hexahedron(id_hex,i, j, k+1);
D3vector cEST = blockgrid->Give_coord_hexahedron(id_hex,i+1,j ,k+1);
D3vector cWNT = blockgrid->Give_coord_hexahedron(id_hex,i, j+1,k+1);
D3vector cENT = blockgrid->Give_coord_hexahedron(id_hex,i+1,j+1,k+1);
// bounding box calculation
D3vector boxWSD, boxENT;
boxWSD.x = MIN(MIN(MIN(cWSD.x,cESD.x),MIN(cWND.x,cEND.x)),
MIN(MIN(cWST.x,cEST.x),MIN(cWNT.x,cENT.x))) - factor *hx;
boxWSD.y = MIN(MIN(MIN(cWSD.y,cESD.y),MIN(cWND.y,cEND.y)),
MIN(MIN(cWST.y,cEST.y),MIN(cWNT.y,cENT.y))) - factor *hy;
boxWSD.z = MIN(MIN(MIN(cWSD.z,cESD.z),MIN(cWND.z,cEND.z)),
MIN(MIN(cWST.z,cEST.z),MIN(cWNT.z,cENT.z))) - factor *hz;
boxENT.x = MAX(MAX(MAX(cWSD.x,cESD.x),MAX(cWND.x,cEND.x)),
MAX(MAX(cWST.x,cEST.x),MAX(cWNT.x,cENT.x))) + factor *hx;
boxENT.y = MAX(MAX(MAX(cWSD.y,cESD.y),MAX(cWND.y,cEND.y)),
MAX(MAX(cWST.y,cEST.y),MAX(cWNT.y,cENT.y))) + factor *hy;
boxENT.z = MAX(MAX(MAX(cWSD.z,cESD.z),MAX(cWND.z,cEND.z)),
MAX(MAX(cWST.z,cEST.z),MAX(cWNT.z,cENT.z))) + factor *hz;
// calculation of indices of a collection of cells of structured grid which contains bounding box
int ilmin = Ganzzahliger_Anteil((boxWSD.x - pWSD.x) / hx);
int jlmin = Ganzzahliger_Anteil((boxWSD.y - pWSD.y) / hy);
int klmin = Ganzzahliger_Anteil((boxWSD.z - pWSD.z) / hz);
int ilmax = Ganzzahliger_Anteil((boxENT.x - pWSD.x) / hx);
int jlmax = Ganzzahliger_Anteil((boxENT.y - pWSD.y) / hy);
int klmax = Ganzzahliger_Anteil((boxENT.z - pWSD.z) / hz);
/*
cout << " indices: "
<< " ilmin: " << ilmin
<< " jlmin: " << jlmin
<< " klmin: " << klmin
<< " ilmax: " << ilmax
<< " jlmax: " << jlmax
<< " klmax: " << klmax
<< " boxWSD.x: " << boxWSD.x
<< " cWSD.x: " << cWSD.x
<< " Nx: " << Nx
<< endl;
*/
/*
bool now;
if(boxWSD.z < 0 && boxENT.z > 0.0 && boxWSD.y < 0.5 && boxENT.y > 0.5 && boxWSD.x < 1.0 && boxENT.x > 1.0 ) {
cout << "\n \n WSD : "; boxWSD.Print();
......@@ -250,7 +242,7 @@ if(boxWSD.z < 0 && boxENT.z > 0.0 && boxWSD.y < 0.5 && boxENT.y > 0.5 && boxWSD.
}
else now = false;
cout << " tt: " << boxWSD.x << " " << pWSD.x << " " << hx << endl;
cout << (boxWSD.x - pWSD.x) << endl;
......@@ -259,101 +251,102 @@ if(boxWSD.z < 0 && boxENT.z > 0.0 && boxWSD.y < 0.5 && boxENT.y > 0.5 && boxWSD.
cout << " z: " << 5.1 << " g: " << Ganzzahliger_Anteil(5.1) << endl;
cout << " z: " << -5.1 << " g: " << Ganzzahliger_Anteil(-5.1) << endl;
*/
if(ilmin<0) ilmin=0;
if(jlmin<0) jlmin=0;
if(klmin<0) klmin=0;
for(int il = ilmin; (il <= ilmax) && (il < nx_);++il)
for(int jl = jlmin; (jl <= jlmax) && (jl < ny_);++jl)
for(int kl = klmin; (kl <= klmax) && (kl < nz_);++kl) {
ploc = D3vector(il * hx, jl * hy, kl * hz) + pWSD;
// cout << "HI" << endl;
typ = -1;
lam = lambda_of_p_in_tet(ploc,cWND,cWNT,cWST,cEST);
if(contained_in_tet(lam)) typ=0;
else {
lam = lambda_of_p_in_tet(ploc,cEST,cWND,cWST,cESD);
if(contained_in_tet(lam)) typ=1;
else {
lam = lambda_of_p_in_tet(ploc,cWND,cWSD,cWST,cESD);
if(contained_in_tet(lam)) typ=2;
else {
lam = lambda_of_p_in_tet(ploc,cEST,cWND,cESD,cEND);
if(contained_in_tet(lam)) typ=3;
else {
lam = lambda_of_p_in_tet(ploc,cENT,cWNT,cEST,cEND);
if(contained_in_tet(lam)) typ=4;
else {
lam = lambda_of_p_in_tet(ploc,cWNT,cWND,cEST,cEND);
if(contained_in_tet(lam)) typ=5;
}
}
}
}
}
/*
cout << " typ " << typ << id_hex
<< " il: " << il
<< " jl: " << jl
<< " kl: " << kl
<< endl;
*/
if(typ!=-1) {
int ind_global;
ind_global = il+nx*(jl+ny*kl);
bool stop;
stop=false;
if(ids_hex[ind_global]!=-1) {
stop=new_lam_worse(lambda[ind_global],lam);
}
if(stop==false) {
ids_hex[ind_global] = id_hex;
ids_i[ind_global] = i;
ids_j[ind_global] = j;
ids_k[ind_global] = k;
typ_tet[ind_global] = typ;
lambda[ind_global] = lam;
}
//go_on = false;
}
*/
if(ilmin<0) ilmin=0;
if(jlmin<0) jlmin=0;
if(klmin<0) klmin=0;
for(int il = ilmin; (il <= ilmax) && (il < nx_);++il)
for(int jl = jlmin; (jl <= jlmax) && (jl < ny_);++jl)
for(int kl = klmin; (kl <= klmax) && (kl < nz_);++kl) {
D3vector ploc = D3vector(il * hx, jl * hy, kl * hz) + pWSD;
// cout << "HI" << endl;
int typ = -1;
D3vector lam = lambda_of_p_in_tet(ploc,cWND,cWNT,cWST,cEST);
if(contained_in_tet(lam)) typ=0;
else {
lam = lambda_of_p_in_tet(ploc,cEST,cWND,cWST,cESD);
if(contained_in_tet(lam)) typ=1;
else {
lam = lambda_of_p_in_tet(ploc,cWND,cWSD,cWST,cESD);
if(contained_in_tet(lam)) typ=2;
else {
lam = lambda_of_p_in_tet(ploc,cEST,cWND,cESD,cEND);
if(contained_in_tet(lam)) typ=3;
else {
lam = lambda_of_p_in_tet(ploc,cENT,cWNT,cEST,cEND);
if(contained_in_tet(lam)) typ=4;
else {
lam = lambda_of_p_in_tet(ploc,cWNT,cWND,cEST,cEND);
if(contained_in_tet(lam)) typ=5;
}
}
}
}
}
/*
cout << " typ " << typ << id_hex
<< " il: " << il
<< " jl: " << jl
<< " kl: " << kl
<< endl;
*/
if(typ!=-1) {
int ind_global;
ind_global = il+nx*(jl+ny*kl);
bool stop;
stop=false;
if(ids_hex[ind_global]!=-1) {
stop=new_lam_worse(lambda[ind_global],lam);
}
/*
cout << " out "
<< " ilmin: " << ilmin
<< " ilmax: " << ilmax
<< " jlmin: " << jlmin
<< " jlmax: " << jlmax
<< " klmin: " << klmin
<< " klmax: " << klmax;
cout << "\n "; cWSD.Print();
cout << "\n "; cESD.Print();
cout << "\n "; cWND.Print();
cout << "\n "; cEND.Print();
cout << "\n "; cWST.Print();
cout << "\n "; cEST.Print();
cout << "\n "; cWNT.Print();
cout << "\n "; cENT.Print();
cout << "\n p: "; ploc.Print();
#pragma omp critical
if(stop==false) {
ids_hex[ind_global] = id_hex;
ids_i[ind_global] = i;
ids_j[ind_global] = j;
ids_k[ind_global] = k;
cout << "\n : "; boxWSD.Print();
cout << "\n : "; boxENT.Print();
typ_tet[ind_global] = typ;
cout << endl;
*/
}
}
lambda[ind_global] = lam;
}
//go_on = false;
}
/*
cout << " out "
<< " ilmin: " << ilmin
<< " ilmax: " << ilmax
<< " jlmin: " << jlmin
<< " jlmax: " << jlmax
<< " klmin: " << klmin
<< " klmax: " << klmax;
cout << "\n "; cWSD.Print();
cout << "\n "; cESD.Print();
cout << "\n "; cWND.Print();
cout << "\n "; cEND.Print();
cout << "\n "; cWST.Print();
cout << "\n "; cEST.Print();
cout << "\n "; cWNT.Print();
cout << "\n "; cENT.Print();
cout << "\n p: "; ploc.Print();
cout << "\n : "; boxWSD.Print();
cout << "\n : "; boxENT.Print();
cout << endl;
*/
}
}
}
......@@ -361,9 +354,9 @@ if(boxWSD.z < 0 && boxENT.z > 0.0 && boxWSD.y < 0.5 && boxENT.y > 0.5 && boxWSD.
if(ids_hex[i]==-1) {
// wir nehmen default value!!
/*
cout << i
<< " Error: Interpolate_on_structured_grid: I cannot interpolate all data!"
<< endl;
cout << i
<< " Error: Interpolate_on_structured_grid: I cannot interpolate all data!"
<< endl;
ids_hex[i] = 0;
*/
}
......@@ -373,6 +366,11 @@ if(boxWSD.z < 0 && boxENT.z > 0.0 && boxWSD.y < 0.5 && boxENT.y > 0.5 && boxWSD.
}
}
/////////////////////////////////////////////////////////////
// 2. Interpolate from blockgrid to blockgrid
/////////////////////////////////////////////////////////////
Interpolate_on_structured_grid::Interpolate_on_structured_grid(int nx_, int ny_, int nz_,
Blockgrid& blockgrid_) {
int Nx, Ny, Nz;
......@@ -665,9 +663,9 @@ Interpolate_on_block_grid::Interpolate_on_block_grid(int nx_, int ny_, int nz_,
data = new double[nx*ny*nz];
hx = (pENT.x - pWSD.x) / (nx+1);
hy = (pENT.y - pWSD.y) / (ny+1); // CHANGES OF RALL MISSING; HAS TO BE FIXED. -> still old version in interpol.cc file used!
hz = (pENT.z - pWSD.z) / (nz+1);
hx = (pENT.x - pWSD.x) / (nx-1);
hy = (pENT.y - pWSD.y) / (ny-1);
hz = (pENT.z - pWSD.z) / (nz-1);
/*
......
program/source/interpol/interpol.h
View file @ b000b0d4
......@@ -68,6 +68,7 @@ class Interpolate_on_structured_grid {
Blockgrid& blockgrid_);
Interpolate_on_structured_grid(int nx_, int ny_, int nz_,
Blockgrid& blockgrid_);
~Interpolate_on_structured_grid();
/***
* interpolates from Variable u(of unstructured blockgrid) to structured data
@param u: Variable on unstructured blockgrid
......
program/source/math_lib/math_lib.h
View file @ b000b0d4
......@@ -24,7 +24,7 @@
#ifndef MATHLIB_H_
#define MATHLIB_H_
#include <limits>
#include "../../../common_source/mathlib/math_lib_common.h"
......@@ -35,6 +35,9 @@
// 2. a simple 3D matrix class and p in hexahedron test
// 3. geometric operators for 3D vectors
//////////////////////////////////////////////////////////////////////
inline bool isNaN(const double pV) {
return (pV != pV) || (fabs(pV) == std::numeric_limits<double>::infinity());
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment