/********************************************************************************** * Copyright 2010 Christoph Pflaum * Department Informatik Lehrstuhl 10 - Systemsimulation * Friedrich-Alexander Universität Erlangen-Nürnberg * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. * You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0 * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. **********************************************************************************/ #include #include #include "../mympi.h" #include "../abbrevi.h" #include "../parameter.h" #include "../math_lib/math_lib.h" #include "../basics/basic.h" #include "elements.h" #include "parti.h" #include "ug.h" #include "blockgrid.h" #define HWSD hex->Give_coord(WSDdir3D) #define HESD hex->Give_coord(ESDdir3D) #define HWND hex->Give_coord(WNDdir3D) #define HEND hex->Give_coord(ENDdir3D) #define HWST hex->Give_coord(WSTdir3D) #define HEST hex->Give_coord(ESTdir3D) #define HWNT hex->Give_coord(WNTdir3D) #define HENT hex->Give_coord(ENTdir3D) int Blockgrid::id_count_grid = 0; Blockgrid::Blockgrid() { id_of_grid = id_count_grid; ++id_count_grid; phase_shift = NULL; bg_coord = NULL; variable_set = false; } Blockgrid::Blockgrid ( Unstructured_grid *ug_ ) { id_of_grid = id_count_grid; ++id_count_grid; phase_shift = NULL; bg_coord = NULL; ug = ug_; number_points = new int[ug->degree_of_freedom() ]; for ( int i=0;idegree_of_freedom();++i ) { number_points[i] = 10; } variable_set = false; } Blockgrid::Blockgrid ( Unstructured_grid *ug_, int N ) { id_of_grid = id_count_grid; ++id_count_grid; phase_shift = NULL; bg_coord = NULL; ug = ug_; number_points = new int[ug->degree_of_freedom() ]; for ( int i=0;idegree_of_freedom();++i ) { number_points[i] = N; } variable_set = false; } Blockgrid::Blockgrid ( Unstructured_grid *ug_, int Nx, int Ny, int Nz ) { id_of_grid = id_count_grid; ++id_count_grid; phase_shift = NULL; bg_coord = NULL; ug = ug_; if ( ug->degree_of_freedom() <3 ) cout << " error in Blockgrid::Blockgrid(Unstructured_grid *ug_, int Nx, int Ny, int Nz) " << endl; number_points = new int[ug->degree_of_freedom() ]; number_points[0] = Nx; number_points[1] = Ny; for ( int i = 2; i < ug->degree_of_freedom();++i ) { number_points[i] = Nz; } variable_set = false; } Blockgrid::Blockgrid ( Unstructured_grid *ug_, int Nx, int Ny, int Nz, int Nr ) { id_of_grid = id_count_grid; ++id_count_grid; phase_shift = NULL; bg_coord = NULL; ug = ug_; if ( ug->degree_of_freedom() !=4 ) { cout << " error in Blockgrid::Blockgrid(Unstructured_grid *ug_,int Nx,int Ny,int Nz,int Nr) " << endl; cout << " number of degree of freedom is: " << ug->degree_of_freedom() << endl; } number_points = new int[4]; number_points[0] = Nx; number_points[1] = Ny; number_points[2] = Nz; number_points[3] = Nr; variable_set = false; } Blockgrid::Blockgrid(Unstructured_grid *ug_, int Nx, int Ny, int Nz, int Nr, int Nr2){ id_of_grid = id_count_grid; ++id_count_grid; phase_shift = NULL; bg_coord = NULL; ug = ug_; if ( ug->degree_of_freedom() !=5 ) { cout << " error in Blockgrid::Blockgrid(Unstructured_grid *ug_,int Nx,int Ny,int Nz,int Nr) " << endl; cout << " number of degree of freedom is: " << ug->degree_of_freedom() << endl; } number_points = new int[5]; number_points[0] = Nx; number_points[1] = Ny; number_points[2] = Nz; number_points[3] = Nr; number_points[4] = Nr2; variable_set = false; } Blockgrid::Blockgrid ( Unstructured_grid *ug_, int Nx, int Ny, std::vector& Nz ) { id_of_grid = id_count_grid; ++id_count_grid; phase_shift = NULL; bg_coord = NULL; ug = ug_; if ( ug->degree_of_freedom() !=Nz.size() +2 ) { std::stringstream ss; ss << " error in Blockgrid::Blockgrid(Unstructured_grid *ug_,int Nx,int Ny,vector& Nz) " << endl << "Nz.size+2 = "<< Nz.size() +2 << endl << "degree_of_freedom = " << ug->degree_of_freedom() << endl; throw ( std::runtime_error ( ss.str() ) ); } number_points = new int[Nz.size() +2]; number_points[0] = Nx; number_points[1] = Ny; // std::cout << "Nx = " << Nx << endl; // std::cout << "Nx = " << Ny << endl; int tmp_num; for ( unsigned int i = 0; i < Nz.size();++i ) { tmp_num = Nz.at ( i ); if ( tmp_num <= 0 ) { std::stringstream ss; ss << " error in Blockgrid::Blockgrid(Unstructured_grid *ug_,int Nx,int Ny,vector& Nz) " << endl << "Nz.at("<< i << ") = " << tmp_num <<" <= 0" << endl; throw ( std::runtime_error ( ss.str() ) ); } // cout << "Nz.at("<< i << ") = " << tmp_num << endl; number_points[i+2] = tmp_num; } variable_set = false; } Blockgrid::Blockgrid ( Unstructured_grid *ug_,std::vector Nvec) { id_of_grid = id_count_grid; ++id_count_grid; phase_shift = NULL; bg_coord = NULL; ug = ug_; if ( ug->degree_of_freedom() !=Nvec.size() ) { std::stringstream ss; ss << " error in Blockgrid::Blockgrid(Unstructured_grid *ug_,std::vector Nvec) " << endl << "Nvec.size() = "<< Nvec.size() << endl << "degree_of_freedom = " << ug->degree_of_freedom() << endl; throw ( std::runtime_error ( ss.str() ) ); } number_points = new int[Nvec.size()]; for(unsigned int i = 0; i < Nvec.size();++i) { number_points[i] = Nvec.at(i); } variable_set = false; } Blockgrid::Blockgrid ( Unstructured_grid *ug_, int Nx, int Ny, std::vector& Nz, int Nr ) { id_of_grid = id_count_grid; ++id_count_grid; phase_shift = NULL; bg_coord = NULL; ug = ug_; if ( ug->degree_of_freedom() !=Nz.size() +3 ) { std::stringstream ss; ss << " error in Blockgrid::Blockgrid(Unstructured_grid *ug_,int Nx,int Ny,vector& Nz,, int Nr) " << endl << "Nz.size+3 = "<< Nz.size() +3 << endl << "degree_of_freedom = " << ug->degree_of_freedom() << endl; throw ( std::runtime_error ( ss.str() ) ); } number_points = new int[Nz.size() +3]; number_points[0] = Nx; number_points[1] = Ny; number_points[Nz.size()+2] = Nr; // std::cout << "Nx = " << Nx << endl; // std::cout << "Nx = " << Ny << endl; int tmp_num; for ( unsigned int i = 0; i < Nz.size();++i ) { tmp_num = Nz.at ( i ); if ( tmp_num <= 0 ) { std::stringstream ss; ss << " error in Blockgrid::Blockgrid(Unstructured_grid *ug_,int Nx,int Ny,vector& Nz,int Nr) " << endl << "Nz.at("<< i << ") = " << tmp_num <<" <= 0" << endl; throw ( std::runtime_error ( ss.str() ) ); } // cout << "Nz.at("<< i << ") = " << tmp_num << endl; number_points[i+2] = tmp_num; } variable_set = false; } Blockgrid::Blockgrid(Unstructured_grid *ug_, std::vector& Nx, std::vector& Ny, std::vector& Nz) { id_of_grid = id_count_grid; ++id_count_grid; phase_shift = NULL; bg_coord = NULL; ug = ug_; int size = Nx.size()+Ny.size()+Nz.size(); if ( ug->degree_of_freedom() != size ) { std::stringstream ss; ss << " error in Blockgrid::Blockgrid(Unstructured_grid *ug_,std::vector& Nx,std::vector& Ny,vector& Nz,, int Nr) " << endl << "Nx.size+Ny.size+Nz.size = "<< Nx.size()+Ny.size()+Nz.size() << endl << "degree_of_freedom = " << ug->degree_of_freedom() << endl; throw ( std::runtime_error ( ss.str() ) ); } number_points = new int[size]; number_points[0] = Nx[0]; number_points[1] = Ny[0]; number_points[2] = Nz[0]; for(int i = 1; i < Nx.size(); ++i) { number_points[i+2] = Nx[i]; } for(int i = 1; i < Ny.size(); ++i) { number_points[Nx.size()+i+1] = Ny[i]; } for(int i = 1; i < Nz.size(); ++i) { number_points[Nx.size()+Ny.size()+i] = Nz[i]; } } Blockgrid::~Blockgrid() { if ( number_points!=NULL ) delete [] number_points; number_points = NULL; if (bg_coord != NULL) delete bg_coord; } int Blockgrid::Give_Nx_hexahedron ( int id ) const { return number_points[ug->Give_edge ( ( ug->Give_hexahedron ( id )-> Give_id_edge ( SDed ) ) )->Give_color_of_edge() ]; } int Blockgrid::Give_Ny_hexahedron ( int id ) const { return number_points[ug->Give_edge ( ( ug->Give_hexahedron ( id )-> Give_id_edge ( WDed ) ) )->Give_color_of_edge() ]; } int Blockgrid::Give_Nz_hexahedron ( int id ) const { return number_points[ug->Give_edge ( ( ug->Give_hexahedron ( id )-> Give_id_edge ( SWed ) ) )->Give_color_of_edge() ]; } int Blockgrid::Give_N_total_hexahedron ( int id ) const { return ( Give_Nx_hexahedron ( id ) +1 ) * ( Give_Ny_hexahedron ( id ) +1 ) * ( Give_Nz_hexahedron ( id ) +1 ); } int Blockgrid::Give_Nx_quadrangle ( int id ) const { return number_points[ug->Give_edge ( ( ug->Give_quadrangle ( id )-> Give_id_edge ( Sdir2D ) ) )->Give_color_of_edge() ]; } int Blockgrid::Give_Ny_quadrangle ( int id ) const { return number_points[ug->Give_edge ( ( ug->Give_quadrangle ( id )-> Give_id_edge ( Wdir2D ) ) )->Give_color_of_edge() ]; } int Blockgrid::Give_N_total_quadrangle ( int id ) const { return ( Give_Nx_quadrangle ( id ) +1 ) * ( Give_Ny_quadrangle ( id ) +1 ); } int Blockgrid::Give_Nx_edge ( int id ) const { return number_points[ug->Give_edge ( id )->Give_color_of_edge() ]; } int Blockgrid::Give_N_color ( int color ) const { if ( developer_version ) if ( color>=ug->degree_of_freedom() || color<0 ) cout << " error in Blockgrid::Give_N_color " << endl; return number_points[color]; } inline double change ( bool ori, double q ) { if ( ori ) return q; return 1.0-q; } inline double find_p ( Edges_cell ed, double eta, double xi,double phi ) { if ( ed<=NTed ) return eta; if ( ed<=ETed ) return xi; if ( developer_version ) if ( ed>NEed ) cout << " error in find_p!" << endl; return phi; } inline std::vector find_p_quad ( dir3D quad, Hexahedron_el *hex ) { // enum Edges_cell { SDed, NDed, STed, NTed, // WDed, EDed, WTed, ETed, // SWed, SEed, NWed, NEed }; std::vector ijk = {1 , 1 , 1}; if (quad == Wdir3D) { hex->orientation[4]; hex->orientation[6]; hex->orientation[8]; hex->orientation[10]; ijk[0] = 0; if (!hex->orientation[4] && !hex->orientation[6]) { ijk[1] = -1; } if (!hex->orientation[8] && !hex->orientation[10]) { ijk[2] = -1; } } else if (quad == Edir3D) { hex->orientation[5]; hex->orientation[7]; hex->orientation[9]; hex->orientation[11]; ijk[0] = 0; if (!hex->orientation[5] && !hex->orientation[7]) { ijk[1] = -1; } if (!hex->orientation[9] && !hex->orientation[11]) { ijk[2] = -1; } } else if (quad == Sdir3D) { hex->orientation[0]; hex->orientation[2]; hex->orientation[8]; hex->orientation[9]; ijk[1] = 0; if (!hex->orientation[0] && !hex->orientation[2]) { ijk[0] = -1; } if (!hex->orientation[8] && !hex->orientation[9]) { ijk[2] = -1; } } else if (quad == Ndir3D) { hex->orientation[1]; hex->orientation[3]; hex->orientation[10]; hex->orientation[11]; ijk[1] = 0; if (!hex->orientation[1] && !hex->orientation[3]) { ijk[0] = -1; } if (!hex->orientation[10] && !hex->orientation[11]) { ijk[2] = -1; } } else if (quad == Ddir3D) { hex->orientation[0]; hex->orientation[1]; hex->orientation[4]; hex->orientation[5]; ijk[2] = 0; if (!hex->orientation[0] && !hex->orientation[1]) { ijk[0] = -1; } if (!hex->orientation[4] && !hex->orientation[5]) { ijk[1] = -1; } } else if (quad == Tdir3D) { hex->orientation[2]; hex->orientation[3]; hex->orientation[6]; hex->orientation[7]; ijk[2] = 0; if (!hex->orientation[2] && !hex->orientation[3]) { ijk[0] = -1; } if (!hex->orientation[6] && !hex->orientation[7]) { ijk[1] = -1; } } return ijk; } D3vector Blockgrid::Give_coord_hexahedron ( int id_hex, int i, int j, int k ) const { if(bg_coord != NULL) if(bg_coord->blockgrid_edge_coordinates_calculated) { int Nx = Give_Nx_hexahedron ( id_hex ) +1; int Ny = Give_Ny_hexahedron ( id_hex ) +1; //cout << "id_hex " << id_hex << " i " << i << " j " << j << " k " << k << endl; //cout << "index " << i + Nx* (j + k * (Ny))<< endl; return bg_coord->blockgrid_hexa_coordinates.at(id_hex).at(i + Nx* (j + k * (Ny)) ); } double eta,xi,phi; D3vector Psi[12]; Hexahedron_el *hex; D3vector PL, PR; D3vector PSW, PSE, PNW, PNE; D3vector Pres, PT, PD; double p_EW, p_NS; double p; eta = ( double ) i/ ( double ) Give_Nx_hexahedron ( id_hex ); xi = ( double ) j/ ( double ) Give_Ny_hexahedron ( id_hex ); phi = ( double ) k/ ( double ) Give_Nz_hexahedron ( id_hex ); hex = ug->Give_hexahedron ( id_hex ); std::vector p_quad_Wdir3D = find_p_quad(Wdir3D,hex); std::vector p_quad_Edir3D = find_p_quad(Edir3D,hex); std::vector p_quad_Sdir3D = find_p_quad(Sdir3D,hex); std::vector p_quad_Ndir3D = find_p_quad(Ndir3D,hex); std::vector p_quad_Ddir3D = find_p_quad(Ddir3D,hex); std::vector p_quad_Tdir3D = find_p_quad(Tdir3D,hex); int idW = hex->Give_id_quadrangle(Wdir3D); int idE = hex->Give_id_quadrangle(Edir3D); int idS = hex->Give_id_quadrangle(Sdir3D); int idN = hex->Give_id_quadrangle(Ndir3D); int idD = hex->Give_id_quadrangle(Ddir3D); int idT = hex->Give_id_quadrangle(Tdir3D); D3vector QS, QN, QE, QW, QT, QD; D3vector QPsi[6]; if(ug->Give_transform_From_Quadrangle()) { //geht bisher nur fuer 1 hexaeder, allgemien so richtig? // idW = 0 ; idE = 1 ; idS = 2 ; idN = 3 ; idD = 4 ; idT = 5; int nEdgeX = Give_Nx_hexahedron ( id_hex ); int nEdgeY = Give_Ny_hexahedron ( id_hex ); int nEdgeZ = Give_Nz_hexahedron ( id_hex ); int quadSaddI = (p_quad_Sdir3D[0] == 1) ? 0 : nEdgeX; int quadSaddK = (p_quad_Sdir3D[2] == 1) ? 0 : nEdgeZ; int quadNaddI = (p_quad_Ndir3D[0] == 1) ? 0 : nEdgeX; int quadNaddK = (p_quad_Ndir3D[2] == 1) ? 0 : nEdgeZ; int quadDaddI = (p_quad_Ddir3D[0] == 1) ? 0 : nEdgeX; int quadDaddJ = (p_quad_Ddir3D[1] == 1) ? 0 : nEdgeY; int quadTaddI = (p_quad_Tdir3D[0] == 1) ? 0 : nEdgeX; int quadTaddJ = (p_quad_Tdir3D[1] == 1) ? 0 : nEdgeY; int quadEaddJ = (p_quad_Edir3D[1] == 1) ? 0 : nEdgeY; int quadEaddK = (p_quad_Edir3D[2] == 1) ? 0 : nEdgeZ; int quadWaddJ = (p_quad_Wdir3D[1] == 1) ? 0 : nEdgeY; int quadWaddK = (p_quad_Wdir3D[2] == 1) ? 0 : nEdgeZ; //with modification QPsi[0] = Give_coord_quadrangle( idW , quadWaddJ + j*p_quad_Wdir3D[1] , quadWaddK + k*p_quad_Wdir3D[2]); QPsi[1] = Give_coord_quadrangle( idE , quadEaddJ + j*p_quad_Edir3D[1] , quadEaddK + k*p_quad_Edir3D[2]); QPsi[2] = Give_coord_quadrangle( idS , quadSaddI + i*p_quad_Sdir3D[0] , quadSaddK + k*p_quad_Sdir3D[2]); QPsi[3] = Give_coord_quadrangle( idN , quadNaddI + i*p_quad_Ndir3D[0] , quadNaddK + k*p_quad_Ndir3D[2]); QPsi[4] = Give_coord_quadrangle( idD , quadDaddI + i*p_quad_Ddir3D[0] , quadDaddJ + j*p_quad_Ddir3D[1]); QPsi[5] = Give_coord_quadrangle( idT , quadTaddI + i*p_quad_Tdir3D[0] , quadTaddJ + j*p_quad_Tdir3D[1]); // Version: 1 x Flaechen + 2 x Kanten + 3 x Ecken // ///////////////////////////////////////////////////////// Pres = (QPsi[0] + eta * (QPsi[1]-QPsi[0]) + QPsi[2] + xi * (QPsi[3]-QPsi[2]) + QPsi[4] + phi * (QPsi[5]-QPsi[4])); //cout << "PresQuad V1 " ; PresQuad.Print();cout<transform[ed] ( change ( hex->orientation[ed],p ) , hex->shiftPointer(ug->Get_pointer_global_data())); PL = hex->Give_coord ( Transform ( ( Edges_cell ) ed,Ldir1D ) ); PR = hex->Give_coord ( Transform ( ( Edges_cell ) ed,Rdir1D ) ); Psi[ed] = Psi[ed] + PL + ( PR-PL ) * p; } Pres = D3vector ( 0.0,0.0,0.0 ); for ( int md = 0; md < 3; ++md ) { if ( md==0 ) // EW { PSW = Psi[SDed]; PSE = Psi[NDed]; PNW = Psi[STed]; PNE = Psi[NTed]; p_EW = xi; p_NS = phi; } if ( md==1 ) // NS { PSW = Psi[WDed]; PSE = Psi[EDed]; PNW = Psi[WTed]; PNE = Psi[ETed]; p_EW = eta; p_NS = phi; } if ( md==2 ) // TD { PSW = Psi[SWed]; PSE = Psi[SEed]; PNW = Psi[NWed]; PNE = Psi[NEed]; p_EW = eta; p_NS = xi; } // Version: 0 x Flaechen + 1 x Kanten + 3 x Ecken // ///////////////////////////////////////////////////////// Pres = Pres + PSW + ( PSE - PSW ) * p_EW + ( PNW - PSW ) * p_NS + ( PNE - PSE - PNW + PSW ) * p_EW * p_NS; } PD = HWSD + ( HESD - HWSD ) * eta + ( HWND - HWSD ) * xi + ( HEND - HESD - HWND + HWSD ) * eta * xi; PT = HWST + ( HEST - HWST ) * eta + ( HWNT - HWST ) * xi + ( HENT - HEST - HWNT + HWST ) * eta * xi; Pres = Pres - 2.0 * ( PD + ( PT-PD ) * phi ); } if ( hex->transform_hex== NULL ) return Pres; else return Pres + hex->transform_hex ( id_hex,eta,xi,phi ); } D3vector Blockgrid::Give_coord_quadrangle ( int id, int i, int j ) const { double eta, xi; D3vector PSW, PSE, PNW, PNE, PW, PE, PN, PS; D3vector PTransFace, PWithoutTrans; Quadrangle_el* quad; Edge_el* edge; quad = this->ug->Give_quadrangle ( id ); eta = ( double ) i/ ( double ) Give_Nx_quadrangle ( id ); xi = ( double ) j/ ( double ) Give_Ny_quadrangle ( id ); // @todo : so sollte es sein PSW = quad->Give_coord ( SWdir2D ); PSE = quad->Give_coord ( SEdir2D ); PNW = quad->Give_coord ( NWdir2D ); PNE = quad->Give_coord ( NEdir2D ); /* // aber ist noch so PSW = ug->Give_point(quad->Give_id_corner(SWdir2D))->Give_coordinate(); PSE = ug->Give_point(quad->Give_id_corner(SEdir2D))->Give_coordinate(); PNW = ug->Give_point(quad->Give_id_corner(NWdir2D))->Give_coordinate(); PNE = ug->Give_point(quad->Give_id_corner(NEdir2D))->Give_coordinate(); */ if ( developer_version ) if ( quad->Give_id_corner ( SWdir2D ) >=quad->Give_id_corner ( SEdir2D ) || quad->Give_id_corner ( SEdir2D ) >=quad->Give_id_corner ( NWdir2D ) || quad->Give_id_corner ( SWdir2D ) >=quad->Give_id_corner ( NEdir2D ) ) cout << " error in Blockgrid::Give_coord_quadrangle! " << endl; //Phillip version with deformed faces if (quad->transform != NULL) { PTransFace = quad->transform ( eta, xi, quad->shiftPointer(ug->Get_pointer_global_data())); PTransFace = PSW + (PSE-PSW) * eta + (PNW-PSW)*xi + (PNE-PSE-PNW+PSW)*xi*eta + PTransFace; return PTransFace; } //Christoph version with deformed edges edge = ug->Give_edge ( quad->Give_id_edge ( Wdir2D ) ); PW = edge->transform ( xi, edge->shiftPointer(ug->Get_pointer_global_data())); PW = PW + PSW + ( PNW-PSW ) * xi; edge = ug->Give_edge ( quad->Give_id_edge ( Sdir2D ) ); PS = edge->transform ( eta, edge->shiftPointer(ug->Get_pointer_global_data())); PS = PS + PSW + ( PSE-PSW ) * eta; edge = ug->Give_edge ( quad->Give_id_edge ( Edir2D ) ); // @todo: warum braucht man hier kein orientation PE = edge->transform ( xi, edge->shiftPointer(ug->Get_pointer_global_data())); PE = PE + PSE + ( PNE-PSE ) * xi; edge = ug->Give_edge ( quad->Give_id_edge ( Ndir2D ) ); // @todo: warum braucht man hier kein orientation PN = edge->transform ( eta , edge->shiftPointer(ug->Get_pointer_global_data())); PN = PN + PNW + ( PNE-PNW ) * eta; return PW + eta* ( PE-PW ) + PS + xi* ( PN-PS ) - ( PSW + ( PSE - PSW ) * eta + ( PNW - PSW ) * xi + ( PNE - PSE - PNW + PSW ) * eta * xi ); } D3vector Blockgrid::Give_coord_edge ( int id, int i ) const { if(bg_coord != NULL) if(bg_coord->blockgrid_edge_coordinates_calculated) { return bg_coord->blockgrid_edge_coordinates.at(id).at(i) ; } Quadrangle_el* quad; Hexahedron_el* hex; D3vector Psi; if (ug->Give_transform_From_Quadrangle()) { int idHex; dir3D dirForQuad; int quadId; if (ug->Give_edge_to_quad_id(id) == -1) { // cout << "slow way, once each, id " << id << endl; //findet hexahedron entsprechend zur kanten ID for (int iterHex = 0 ; iterHex < ug->Give_number_hexahedra() ; iterHex++) { idHex = iterHex; hex = ug->Give_hexahedron(iterHex); for (int iter = 0 ; iter < 12 ; iter++) { int idCorner = hex->Give_id_edge(Edges_cell(iter)); if (idCorner == id) { iterHex = ug->Give_number_hexahedra() ; iter = 12; } } } //findet quad ID sowie ausrichtung entsprechend zur kanten ID for ( int iter = 0 ; iter < 6 ; iter++) { quad = ug->Give_quadrangle(hex->Give_id_quadrangle((dir3D)iter)); dirForQuad = (dir3D)iter ; for (int iterEdge = 0 ; iterEdge < 4 ; iterEdge++) { quadId = quad->Give_id_edge((dir2D)iterEdge); if (id == quadId) { quadId = hex->Give_id_quadrangle((dir3D)iter); iterEdge = 4; iter = 6; } } } ug->Set_edge_to_quad_id(id,quadId); ug->Set_edge_to_quad_dir(id,(int)dirForQuad); } else { quadId = ug->Give_edge_to_quad_id(id); dirForQuad = (dir3D)ug->Give_edge_to_quad_dir(id); quad = ug->Give_quadrangle(quadId); } int Nx = Give_Nx_quadrangle(dirForQuad); int Ny = Give_Ny_quadrangle(dirForQuad); if (id == quad->Give_id_edge(Wdir2D)) //Wdir2D { Psi = Give_coord_quadrangle(quadId,0,i); } if (id == quad->Give_id_edge(Edir2D)) //Edir2D { Psi = Give_coord_quadrangle(quadId,Nx,i); } if (id == quad->Give_id_edge(Sdir2D)) //Sdir2D { Psi = Give_coord_quadrangle(quadId,i,0); } if (id == quad->Give_id_edge(Ndir2D)) //Ndir2D { Psi = Give_coord_quadrangle(quadId,i,Ny); } return Psi; } //Christoph version, if quadrangles are NOT defined. double eta; D3vector PL, PR; Edge_el* edge; edge = ug->Give_edge ( id ); eta = ( double ) i/ ( double ) Give_Nx_edge ( id ); //cout << "psi new "; Psi.Print();cout << endl; Psi = edge->transform ( eta , edge->shiftPointer(ug->Get_pointer_global_data())); PL = edge->Give_coord ( Ldir1D ); PR = edge->Give_coord ( Rdir1D ); /* PL = ug->Give_point(edge->Give_id_corner_W())->Give_coordinate(); PR = ug->Give_point(edge->Give_id_corner_E())->Give_coordinate(); */ // cout << "psi old " ; (Psi + PL + ( PR-PL ) * eta).Print();cout << endl; return Psi + PL + ( PR-PL ) * eta; } D3vector Blockgrid::Give_coord_point ( int id_point ) const { return ug->Give_point ( id_point )->Give_coordinate(); } int Blockgrid::Total_number_of_points() const { using namespace std; int num_points, num_points_hex, num_points_quad, num_points_edges; int Nx, Ny, Nz; num_points = num_points_hex = num_points_quad = num_points_edges = 0; // hexahedra for ( int id=0;idGive_number_hexahedra();++id ) { Nx = Give_Nx_hexahedron ( id ); Ny = Give_Ny_hexahedron ( id ); Nz = Give_Nz_hexahedron ( id ); num_points_hex += ( Nx-1 ) * ( Ny-1 ) * ( Nz-1 ); } // quadrangles for ( int id=0;idGive_number_quadrangles();++id ) { Nx = Give_Nx_quadrangle ( id ); Ny = Give_Ny_quadrangle ( id ); num_points_quad += ( Nx-1 ) * ( Ny-1 ); } // edges for ( int id=0;idGive_number_edges();++id ) { Nx = Give_Nx_edge ( id ); num_points_edges += ( Nx-1 ); } // points num_points = ug->Give_number_points(); // alltogether num_points = num_points + num_points_hex + num_points_quad + num_points_edges; /* cout << "total number of points of blockgrid: " << num_points << endl; getchar();*/ return num_points; } void Blockgrid_coordinates::init_blockgrid_coordinates() { blockgrid_edge_coordinates_calculated = false; blockgrid_hexa_coordinates_calculated = false; blockgrid_hexa_boundaries_calculated = false; blockgrid_hexa_coordinates.resize(bg->Give_unstructured_grid()->Give_number_hexahedra()); for (int id_hex = 0 ; id_hex < bg->Give_unstructured_grid()->Give_number_hexahedra() ; id_hex++) { int Nx = bg->Give_Nx_hexahedron(id_hex)+1; int Ny = bg->Give_Ny_hexahedron(id_hex)+1; int Nz = bg->Give_Nz_hexahedron(id_hex)+1; blockgrid_hexa_coordinates.at(id_hex).resize(Nx*Ny*Nz); } for (int id_hex = 0 ; id_hex < bg->Give_unstructured_grid()->Give_number_hexahedra() ; id_hex++) { int Nx = bg->Give_Nx_hexahedron(id_hex)+1; int Ny = bg->Give_Ny_hexahedron(id_hex)+1; int Nz = bg->Give_Nz_hexahedron(id_hex)+1; for(int i=0;iGive_coord_hexahedron(id_hex,i,j,k); } } blockgrid_edge_coordinates.resize(bg->Give_unstructured_grid()->Give_number_edges()); for (int id_edge = 0 ; id_edge < bg->Give_unstructured_grid()->Give_number_edges() ; id_edge++) { int Nx = bg->Give_Nx_edge(id_edge)+1; blockgrid_edge_coordinates.at(id_edge).resize(Nx); } for (int id_edge = 0 ; id_edge < bg->Give_unstructured_grid()->Give_number_edges() ; id_edge++) { int Nx = bg->Give_Nx_edge(id_edge)+1; for(int i=0;iGive_coord_edge(id_edge,i); } } blockgrid_edge_coordinates_calculated = true; blockgrid_hexa_coordinates_calculated = true; } void Blockgrid_coordinates::init_blockgrid_coordinates_boundary() { blockgrid_hexa_boundary.clear(); int counter(0); blockgrid_hexa_boundary.resize(bg->Give_unstructured_grid()->Give_number_hexahedra()); //outer loop for (int id_hex = 0 ; id_hex < bg->Give_unstructured_grid()->Give_number_hexahedra() ; id_hex++) { int Nx = bg->Give_Nx_hexahedron(id_hex)+1; int Ny = bg->Give_Ny_hexahedron(id_hex)+1; int Nz = bg->Give_Nz_hexahedron(id_hex)+1; blockgrid_hexa_boundary.at(id_hex).resize(Nx*Ny*Nz); //for(int i=0;i 0 || i < Nx-1) ^ (j > 0 || j < Ny-1) ^ (k > 0 || k < Nz-1)) //remove true here? if (false || (i == 0 || i == Nx-1) || (j== 0 || j == Ny-1) || (k == 0 || k == Nz-1)) { //inner loop for (int id_hexInner = 0 ; id_hexInner < bg->Give_unstructured_grid()->Give_number_hexahedra() ; id_hexInner++) { if (id_hex == id_hexInner) { //id_hexInner++; } else { //cout << "id_hexInner " << id_hexInner<< endl; int NxInner = bg->Give_Nx_hexahedron(id_hexInner)+1; int NyInner = bg->Give_Ny_hexahedron(id_hexInner)+1; int NzInner = bg->Give_Nz_hexahedron(id_hexInner)+1; //for(int iInner=0;iInner 0 || iInner < NxInner-1) ^ (jInner > 0 || jInner < NyInner-1) ^ (kInner > 0 || kInner < NzInner-1)) if (false || (iInner == 0 || iInner == NxInner-1) || (jInner == 0 || jInner == NyInner-1) || (kInner == 0 || kInner == NzInner-1)) { if(D3VectorNormSquared(bg->Give_coord_hexahedron(id_hex,i,j,k) - bg->Give_coord_hexahedron(id_hexInner,iInner,jInner,kInner)) < 1e-10 ) { int convertedLocalIndex = i +(Nx)*(j +(Ny)* k); counter++; blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(id_hexInner); blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(iInner); blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(jInner); blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(kInner); } } } } } } } } blockgrid_hexa_boundaries_calculated = true; } void Blockgrid_coordinates::init_blockgrid_coordinates_boundary_fast() { blockgrid_hexa_boundary.clear(); counter = 0; blockgrid_hexa_boundary.resize(bg->Give_unstructured_grid()->Give_number_hexahedra()); //outer loop for (int id_hex = 0 ; id_hex < bg->Give_unstructured_grid()->Give_number_hexahedra() ; id_hex++) { int Nx = bg->Give_Nx_hexahedron(id_hex); int Ny = bg->Give_Ny_hexahedron(id_hex); int Nz = bg->Give_Nz_hexahedron(id_hex); blockgrid_hexa_boundary.at(id_hex).resize((Nx+1)*(Ny+1)*(Nz+1)); for(int k=0;k<=Nz;k+=Nz/2) for(int j=0;j<=Ny;j+=Ny/2) for(int i=0;i<=Nx;i+=Nx/2) { bool A = (i == 0 || i == Nx); bool B = (j == 0 || j == Ny); bool C = (k == 0 || k == Nz); if ((A&!B&!C)|(!A&B&!C)|(!A&!B&C)) { //inner loop for (int id_hexInner = 0 ; id_hexInner < bg->Give_unstructured_grid()->Give_number_hexahedra() ; id_hexInner++) { if (id_hex != id_hexInner) { //cout << "id_hexInner " << id_hexInner<< endl; int NxInner = bg->Give_Nx_hexahedron(id_hexInner); int NyInner = bg->Give_Ny_hexahedron(id_hexInner); int NzInner = bg->Give_Nz_hexahedron(id_hexInner); for(int kInner=0;kInner<=NzInner;kInner+=NzInner/2) for(int jInner=0;jInner<=NyInner;jInner+=NyInner/2) for(int iInner=0;iInner<=NxInner;iInner+=NxInner/2) { bool AI = (iInner == 0 || iInner == NxInner); bool BI = (jInner == 0 || jInner == NyInner); bool CI = (kInner == 0 || kInner == NzInner); if ((AI&!BI&!CI)|(!AI&BI&!CI)|(!AI&!BI&CI)) { if(D3VectorNormSquared(bg->Give_coord_hexahedron(id_hex,i,j,k) - bg->Give_coord_hexahedron(id_hexInner,iInner,jInner,kInner)) < 1e-10 ) { //add whole plane : addBoundaryPlane(id_hex,i,j,k,id_hexInner,iInner,jInner,kInner, Nx, Ny, Nz, NxInner, NyInner, NzInner); // int convertedLocalIndex = i +(Nx+1)*(j +(Ny+1)* k); // counter++; // blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(id_hexInner); // blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(iInner); // blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(jInner); // blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(kInner); } } } } } } } } blockgrid_hexa_boundaries_calculated = true; } void Blockgrid_coordinates::addBoundaryPlane(int id_hex, int i, int j, int k, int id_hexInner, int iInner, int jInner, int kInner, int Nx, int Ny, int Nz, int NxInner, int NyInner, int NzInner) { int PlaneHits{0}; int tempTries{0}; if (i == 0 || i == Nx) { //vary j,k for (j = 0 ; j <= Ny;j++)for (k = 0 ; k <= Nz;k++) { if (iInner == 0 || iInner == NxInner) { //vary j,k for (int jInnerT = 0 ; jInnerT <= NyInner;jInnerT++)for (int kInnerT = 0 ; kInnerT <= NzInner;kInnerT++) { tempTries++; if(D3VectorNormSquared(bg->Give_coord_hexahedron(id_hex,i,j,k) - bg->Give_coord_hexahedron(id_hexInner,iInner,jInnerT,kInnerT)) < 1e-10 ) { int convertedLocalIndex = i +(Nx+1)*(j +(Ny+1)* k); counter++; PlaneHits++; blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(id_hexInner); blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(iInner); blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(jInnerT); blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(kInnerT); } } } else if (jInner == 0 || jInner == NyInner) { //vary i,k for (int iInnerT = 0 ; iInnerT <= NxInner;iInnerT++)for (int kInnerT = 0 ; kInnerT <= NzInner;kInnerT++) { tempTries++; if(D3VectorNormSquared(bg->Give_coord_hexahedron(id_hex,i,j,k) - bg->Give_coord_hexahedron(id_hexInner,iInnerT,jInner,kInnerT)) < 1e-10 ) { int convertedLocalIndex = i +(Nx+1)*(j +(Ny+1)* k); counter++; PlaneHits++; blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(id_hexInner); blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(iInnerT); blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(jInner); blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(kInnerT); } } } else { //vary i,j for (int iInnerT = 0 ; iInnerT <= NxInner;iInnerT++)for (int jInnerT = 0 ; jInnerT <= NyInner;jInnerT++) { if(D3VectorNormSquared(bg->Give_coord_hexahedron(id_hex,i,j,k) - bg->Give_coord_hexahedron(id_hexInner,iInnerT,jInnerT,kInner)) < 1e-10 ) { tempTries++; int convertedLocalIndex = i +(Nx+1)*(j +(Ny+1)* k); counter++; PlaneHits++; blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(id_hexInner); blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(iInnerT); blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(jInnerT); blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(kInner); } } } } } else if (j == 0 || j == Ny) { //vary i,k for (i = 0 ; i <= Nx;i++)for (k = 0 ; k <= Nz;k++) { if (iInner == 0 || iInner == NxInner) { //vary j,k for (int jInnerT = 0 ; jInnerT <= NyInner;jInnerT++)for (int kInnerT = 0 ; kInnerT <= NzInner;kInnerT++) { tempTries++; if(D3VectorNormSquared(bg->Give_coord_hexahedron(id_hex,i,j,k) - bg->Give_coord_hexahedron(id_hexInner,iInner,jInnerT,kInnerT)) < 1e-10 ) { int convertedLocalIndex = i +(Nx+1)*(j +(Ny+1)* k); counter++; PlaneHits++; blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(id_hexInner); blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(iInner); blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(jInnerT); blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(kInnerT); } } } else if (jInner == 0 || jInner == NyInner) { //vary i,k for (int iInnerT = 0 ; iInnerT <= NxInner;iInnerT++)for (int kInnerT = 0 ; kInnerT <= NzInner;kInnerT++) { tempTries++; if(D3VectorNormSquared(bg->Give_coord_hexahedron(id_hex,i,j,k) - bg->Give_coord_hexahedron(id_hexInner,iInnerT,jInner,kInnerT)) < 1e-10 ) { int convertedLocalIndex = i +(Nx+1)*(j +(Ny+1)* k); counter++; PlaneHits++; blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(id_hexInner); blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(iInnerT); blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(jInner); blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(kInnerT); } } } else { //vary i,j for (int iInnerT = 0 ; iInnerT <= NxInner;iInnerT++)for (int jInnerT = 0 ; jInnerT <= NyInner;jInnerT++) { tempTries++; if(D3VectorNormSquared(bg->Give_coord_hexahedron(id_hex,i,j,k) - bg->Give_coord_hexahedron(id_hexInner,iInnerT,jInnerT,kInner)) < 1e-10 ) { int convertedLocalIndex = i +(Nx+1)*(j +(Ny+1)* k); counter++; PlaneHits++; blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(id_hexInner); blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(iInnerT); blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(jInnerT); blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(kInner); } } } } } else { //vary i,j for (i = 0 ; i <= Nx;i++)for (j = 0 ; j <= Ny;j++) { if (iInner == 0 || iInner == NxInner) { //vary j,k for (int jInnerT = 0 ; jInnerT <= NyInner;jInnerT++)for (int kInnerT = 0 ; kInnerT <= NzInner;kInnerT++) { tempTries++; if(D3VectorNormSquared(bg->Give_coord_hexahedron(id_hex,i,j,k) - bg->Give_coord_hexahedron(id_hexInner,iInner,jInnerT,kInnerT)) < 1e-10 ) { int convertedLocalIndex = i +(Nx+1)*(j +(Ny+1)* k); counter++; PlaneHits++; blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(id_hexInner); blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(iInner); blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(jInnerT); blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(kInnerT); } } } else if (jInner == 0 || jInner == NyInner) { //vary i,k for (int iInnerT = 0 ; iInnerT <= NxInner;iInnerT++)for (int kInnerT = 0 ; kInnerT <= NzInner;kInnerT++) { tempTries++; if(D3VectorNormSquared(bg->Give_coord_hexahedron(id_hex,i,j,k) - bg->Give_coord_hexahedron(id_hexInner,iInnerT,jInner,kInnerT)) < 1e-10 ) { int convertedLocalIndex = i +(Nx+1)*(j +(Ny+1)* k); counter++; PlaneHits++; blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(id_hexInner); blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(iInnerT); blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(jInner); blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(kInnerT); } } } else { //vary i,j for (int iInnerT = 0 ; iInnerT <= NxInner;iInnerT++)for (int jInnerT = 0 ; jInnerT <= NyInner;jInnerT++) { tempTries++; if(D3VectorNormSquared(bg->Give_coord_hexahedron(id_hex,i,j,k) - bg->Give_coord_hexahedron(id_hexInner,iInnerT,jInnerT,kInner)) < 1e-10 ) { int convertedLocalIndex = i +(Nx+1)*(j +(Ny+1)* k); counter++; PlaneHits++; blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(id_hexInner); blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(iInnerT); blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(jInnerT); blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(kInner); } } } } } // std::cout << "tempTries " << tempTries << std::endl; // std::cout << "PlaneHits " << PlaneHits << std::endl; // std::cout << "PlaneHits " << PlaneHits << std::endl; } void Blockgrid_coordinates::delete_blockgrid_coordinates() { blockgrid_edge_coordinates_calculated = false; blockgrid_hexa_coordinates_calculated = false; blockgrid_hexa_coordinates.clear(); blockgrid_hexa_boundary.clear(); }