Commit 7e6f42ac authored by Christoph Pflaum's avatar Christoph Pflaum
Browse files

blockgrid coordinates in constructor und diffop local stiffness matrices...

blockgrid coordinates in constructor und diffop local stiffness matrices effizienter falls constant in z-Richtung
parent 4f7c2a61
......@@ -41,8 +41,25 @@
int Blockgrid::id_count_grid = 0;
bool Blockgrid::constructBlockCoordinates = false;
void Blockgrid::constructBlockCoordinatesIfRequired() {
if(constructBlockCoordinates) {
Blockgrid_coordinates *bg_coord__ = new Blockgrid_coordinates(this);
Set_blockgrid_coordinates( bg_coord__ );
bg_coord->init_blockgrid_coordinates();
}
else {
Set_blockgrid_coordinates(NULL);
}
}
void Blockgrid::Set_blockgrid_coordinates(Blockgrid_coordinates* bg_coord_ ) {
if(bg_coord == NULL) delete bg_coord;
bg_coord = bg_coord_;
}
Blockgrid::Blockgrid() {
id_of_grid = id_count_grid; ++id_count_grid;
......@@ -56,7 +73,6 @@ Blockgrid::Blockgrid() {
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;i<ug->degree_of_freedom();++i )
......@@ -64,12 +80,13 @@ Blockgrid::Blockgrid ( Unstructured_grid *ug_ ) {
number_points[i] = 10;
}
variable_set = false;
}
bg_coord = NULL;
constructBlockCoordinatesIfRequired();
}
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;i<ug->degree_of_freedom();++i )
......@@ -77,12 +94,13 @@ Blockgrid::Blockgrid ( Unstructured_grid *ug_, int N ) {
number_points[i] = N;
}
variable_set = false;
bg_coord = NULL;
constructBlockCoordinatesIfRequired();
}
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) "
......@@ -95,15 +113,15 @@ Blockgrid::Blockgrid ( Unstructured_grid *ug_, int Nx, int Ny, int Nz ) {
{
number_points[i] = Nz;
}
variable_set = false;
bg_coord = NULL;
constructBlockCoordinatesIfRequired();
}
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) "
......@@ -118,12 +136,13 @@ Blockgrid::Blockgrid ( Unstructured_grid *ug_, int Nx, int Ny, int Nz, int Nr )
number_points[3] = Nr;
variable_set = false;
bg_coord = NULL;
constructBlockCoordinatesIfRequired();
}
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) "
......@@ -139,12 +158,13 @@ Blockgrid::Blockgrid(Unstructured_grid *ug_, int Nx, int Ny, int Nz, int Nr, int
number_points[4] = Nr2;
variable_set = false;
bg_coord = NULL;
constructBlockCoordinatesIfRequired();
}
Blockgrid::Blockgrid ( Unstructured_grid *ug_, int Nx, int Ny, std::vector<int>& 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 )
{
......@@ -178,6 +198,8 @@ Blockgrid::Blockgrid ( Unstructured_grid *ug_, int Nx, int Ny, std::vector<int>&
}
variable_set = false;
bg_coord = NULL;
constructBlockCoordinatesIfRequired();
}
......@@ -186,7 +208,6 @@ Blockgrid::Blockgrid ( Unstructured_grid *ug_, int Nx, int Ny, std::vector<int>&
Blockgrid::Blockgrid ( Unstructured_grid *ug_,std::vector<int> 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() )
{
......@@ -202,6 +223,8 @@ Blockgrid::Blockgrid ( Unstructured_grid *ug_,std::vector<int> Nvec) {
}
variable_set = false;
bg_coord = NULL;
constructBlockCoordinatesIfRequired();
}
......@@ -212,7 +235,6 @@ Blockgrid::Blockgrid ( Unstructured_grid *ug_, int Nx, int Ny, std::vector<int>&
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 )
{
......@@ -246,6 +268,8 @@ Blockgrid::Blockgrid ( Unstructured_grid *ug_, int Nx, int Ny, std::vector<int>&
}
variable_set = false;
bg_coord = NULL;
constructBlockCoordinatesIfRequired();
}
......@@ -254,7 +278,6 @@ Blockgrid::Blockgrid(Unstructured_grid *ug_, std::vector<int>& Nx, std::vector<i
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();
......@@ -283,6 +306,8 @@ Blockgrid::Blockgrid(Unstructured_grid *ug_, std::vector<int>& Nx, std::vector<i
{
number_points[Nx.size()+Ny.size()+i] = Nz[i];
}
bg_coord = NULL;
constructBlockCoordinatesIfRequired();
}
Blockgrid::~Blockgrid()
......@@ -478,10 +503,7 @@ inline std::vector<int> find_p_quad ( dir3D quad, Hexahedron_el *hex )
}
D3vector Blockgrid::Give_coord_hexahedron ( int id_hex,
int i, int j, int k ) const
{
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;
......
......@@ -98,6 +98,8 @@ class Blockgrid {
Blockgrid(Unstructured_grid *ug_, std::vector<int>& Nx, std::vector<int>& Ny, std::vector<int>& Nz);
Blockgrid(Unstructured_grid *ug_, std::vector<int> Nvec);
static bool constructBlockCoordinates;
/** Constructor for vector of Cylinders
*** @param Nx: Number of grid points in x-direction
*** @param Ny: Number of grid points in y-direction
......@@ -111,7 +113,7 @@ class Blockgrid {
void Set_grid_points(int d, int Nd);
Unstructured_grid* Give_unstructured_grid() const { return ug; };
Blockgrid_coordinates* Give_blockgrid_coordinates() const { return bg_coord; };
void Set_blockgrid_coordinates(Blockgrid_coordinates* bg_coord_ ) { bg_coord = bg_coord_; };
void Set_blockgrid_coordinates(Blockgrid_coordinates* bg_coord_ );
D3vector Give_coord_hexahedron(int id, int i, int j, int k) const ;
D3vector Give_coord_quadrangle(int id, int i, int j, bool invert = false) const;
......@@ -169,10 +171,8 @@ class Blockgrid {
private:
Variable<std::complex<double> > *phase_shift;
void constructBlockCoordinatesIfRequired();
bool variable_set; // set true falls eine Variable konstruiert
int id_of_grid;
......
......@@ -140,7 +140,8 @@ class Unstructured_grid : public Partitioning {
public:
Unstructured_grid();
~Unstructured_grid();
TypeOfUG getUgTyp() { return typFuerSlice; }
TypeOfUG getUgTyp() { return typFuerSlice; }
bool isConstantCellInZ() { return (typFuerSlice>0) && (typFuerSlice<=5); }
std::vector<double> getConstructionParameters() { return constructionParameters; }
std::vector<double> getConstructionParametersLx() { return construction_lx; }
......
......@@ -736,94 +736,142 @@ void Local_stiffness_matrix<TYPE2>::Calculate ( const B& bilinear_form ) {
Unstructured_grid* ug = blockgrid->Give_unstructured_grid();
//std::vector<std::vector<TYPE2> > loc_stiff;
//std::vector<double> points ( 24, 0. ); //Länge 24
//_COLSAMM_::ELEMENTS::_Hexahedron_2<_COLSAMM_::Gauss2, TYPE2> G;
for ( id = 0; id < blockgrid->Give_unstructured_grid()->Give_number_hexahedra(); ++id ) {
if ( ug->Give_hexahedron ( id )->my_object ( my_rank ) ) {
for(id = 0; id < blockgrid->Give_unstructured_grid()->Give_number_hexahedra(); ++id ) {
if(ug->Give_hexahedron ( id )->my_object ( my_rank ) ) {
Nx = blockgrid->Give_Nx_hexahedron ( id );
Ny = blockgrid->Give_Ny_hexahedron ( id );
Nz = blockgrid->Give_Nz_hexahedron ( id );
N_total = Nx * Ny * Nz * 64;
Nx = blockgrid->Give_Nx_hexahedron ( id );
Ny = blockgrid->Give_Ny_hexahedron ( id );
Nz = blockgrid->Give_Nz_hexahedron ( id );
N_total = Nx * Ny * Nz * 64;
if(ug->isConstantCellInZ()) {
int k = 0;
#pragma omp parallel for num_threads(UGBlocks::numThreadsToTake) if(UGBlocks::useOpenMP)
for(int j = 0;j < Ny;++j) {
std::vector<std::vector<TYPE2> > loc_stiff;
std::vector<double> points ( 24, 0. ); //Länge 24
_COLSAMM_::ELEMENTS::_Hexahedron_2<_COLSAMM_::Gauss2, TYPE2> G;
//neu eingebaut!!!
#pragma omp parallel for num_threads(UGBlocks::numThreadsToTake) if(UGBlocks::useOpenMP)
for(int k = 0;k < Nz;++k ) {
std::vector<std::vector<TYPE2> > loc_stiff;
std::vector<double> points ( 24, 0. ); //Länge 24
_COLSAMM_::ELEMENTS::_Hexahedron_2<_COLSAMM_::Gauss2, TYPE2> G;
for(int j = 0;j < Ny;++j )
for(int i = 0;i < Nx;++i ) {
if ( developer_version )
if ( Ind_loc_matrix_hexahedra ( i, j, k ) *64 >= N_total )
std::cout << " error in Local_stiffness_matrix::Calculate" << std::endl;
for(int i = 0;i < Nx;++i ) {
if(developer_version) if(Ind_loc_matrix_hexahedra ( i, j, k ) *64 >= N_total )
std::cout << " error in Local_stiffness_matrix::Calculate" << std::endl;
points[0] = blockgrid->Give_coord_hexahedron ( id, i, j, k ).x;
points[1] = blockgrid->Give_coord_hexahedron ( id, i, j, k ).y;
points[2] = blockgrid->Give_coord_hexahedron ( id, i, j, k ).z;
points[3] = blockgrid->Give_coord_hexahedron ( id, i + 1, j, k ).x;
points[4] = blockgrid->Give_coord_hexahedron ( id, i + 1, j, k ).y;
points[5] = blockgrid->Give_coord_hexahedron ( id, i + 1, j, k ).z;
points[6] = blockgrid->Give_coord_hexahedron ( id, i, j + 1, k ).x;
points[7] = blockgrid->Give_coord_hexahedron ( id, i, j + 1, k ).y;
points[0] = blockgrid->Give_coord_hexahedron ( id, i, j, k ).x;
points[1] = blockgrid->Give_coord_hexahedron ( id, i, j, k ).y;
points[2] = blockgrid->Give_coord_hexahedron ( id, i, j, k ).z;
points[8] = blockgrid->Give_coord_hexahedron ( id, i, j + 1, k ).z;
points[3] = blockgrid->Give_coord_hexahedron ( id, i + 1, j, k ).x;
points[4] = blockgrid->Give_coord_hexahedron ( id, i + 1, j, k ).y;
points[5] = blockgrid->Give_coord_hexahedron ( id, i + 1, j, k ).z;
points[9] = blockgrid->Give_coord_hexahedron ( id, i + 1, j + 1, k ).x;
points[6] = blockgrid->Give_coord_hexahedron ( id, i, j + 1, k ).x;
points[7] = blockgrid->Give_coord_hexahedron ( id, i, j + 1, k ).y;
points[8] = blockgrid->Give_coord_hexahedron ( id, i, j + 1, k ).z;
points[10] = blockgrid->Give_coord_hexahedron ( id, i + 1, j + 1, k ).y;
points[9] = blockgrid->Give_coord_hexahedron ( id, i + 1, j + 1, k ).x;
points[10] = blockgrid->Give_coord_hexahedron ( id, i + 1, j + 1, k ).y;
points[11] = blockgrid->Give_coord_hexahedron ( id, i + 1, j + 1, k ).z;
points[11] = blockgrid->Give_coord_hexahedron ( id, i + 1, j + 1, k ).z;
points[12] = blockgrid->Give_coord_hexahedron ( id, i, j, k + 1 ).x;
points[13] = blockgrid->Give_coord_hexahedron ( id, i, j, k + 1 ).y;
points[14] = blockgrid->Give_coord_hexahedron ( id, i, j, k + 1 ).z;
points[12] = blockgrid->Give_coord_hexahedron ( id, i, j, k + 1 ).x;
points[15] = blockgrid->Give_coord_hexahedron ( id, i + 1, j, k + 1 ).x;
points[16] = blockgrid->Give_coord_hexahedron ( id, i + 1, j, k + 1 ).y;
points[17] = blockgrid->Give_coord_hexahedron ( id, i + 1, j, k + 1 ).z;
points[13] = blockgrid->Give_coord_hexahedron ( id, i, j, k + 1 ).y;
points[18] = blockgrid->Give_coord_hexahedron ( id, i, j + 1, k + 1 ).x;
points[19] = blockgrid->Give_coord_hexahedron ( id, i, j + 1, k + 1 ).y;
points[20] = blockgrid->Give_coord_hexahedron ( id, i, j + 1, k + 1 ).z;
points[14] = blockgrid->Give_coord_hexahedron ( id, i, j, k + 1 ).z;
points[21] = blockgrid->Give_coord_hexahedron ( id, i + 1, j + 1, k + 1 ).x;
points[22] = blockgrid->Give_coord_hexahedron ( id, i + 1, j + 1, k + 1 ).y;
points[23] = blockgrid->Give_coord_hexahedron ( id, i + 1, j + 1, k + 1 ).z;
points[15] = blockgrid->Give_coord_hexahedron ( id, i + 1, j, k + 1 ).x;
G(points);
points[16] = blockgrid->Give_coord_hexahedron ( id, i + 1, j, k + 1 ).y;
loc_stiff = G.integrate ( bilinear_form );
points[17] = blockgrid->Give_coord_hexahedron ( id, i + 1, j, k + 1 ).z;
for(int v = 0;v < 8;++v )
for(int u = 0;u < 8;++u )
loc_m_hexahedra[id][Ind_loc_matrix_hexahedra ( i,j,k ) *64+8*v+u]
= loc_stiff[u][v];
}
}
#pragma omp parallel for num_threads(UGBlocks::numThreadsToTake) if(UGBlocks::useOpenMP)
for(int k = 1;k < Nz;++k ) {
for(int j = 0;j < Ny;++j) {
for(int i = 0;i < Nx;++i ) {
for(int vu = 0;vu < 64;++vu )
loc_m_hexahedra[id][Ind_loc_matrix_hexahedra(i,j,k) *64+vu] =
loc_m_hexahedra[id][Ind_loc_matrix_hexahedra ( i,j,0 ) *64+vu];
}
}
}
}
else {
#pragma omp parallel for num_threads(UGBlocks::numThreadsToTake) if(UGBlocks::useOpenMP)
for(int k = 0;k < Nz;++k ) {
std::vector<std::vector<TYPE2> > loc_stiff;
std::vector<double> points ( 24, 0. ); //Länge 24
_COLSAMM_::ELEMENTS::_Hexahedron_2<_COLSAMM_::Gauss2, TYPE2> G;
for(int j = 0;j < Ny;++j )
for(int i = 0;i < Nx;++i ) {
if(developer_version) if(Ind_loc_matrix_hexahedra ( i, j, k ) *64 >= N_total )
std::cout << " error in Local_stiffness_matrix::Calculate" << std::endl;
points[0] = blockgrid->Give_coord_hexahedron ( id, i, j, k ).x;
points[1] = blockgrid->Give_coord_hexahedron ( id, i, j, k ).y;
points[2] = blockgrid->Give_coord_hexahedron ( id, i, j, k ).z;
points[18] = blockgrid->Give_coord_hexahedron ( id, i, j + 1, k + 1 ).x;
points[3] = blockgrid->Give_coord_hexahedron ( id, i + 1, j, k ).x;
points[4] = blockgrid->Give_coord_hexahedron ( id, i + 1, j, k ).y;
points[5] = blockgrid->Give_coord_hexahedron ( id, i + 1, j, k ).z;
points[19] = blockgrid->Give_coord_hexahedron ( id, i, j + 1, k + 1 ).y;
points[6] = blockgrid->Give_coord_hexahedron ( id, i, j + 1, k ).x;
points[7] = blockgrid->Give_coord_hexahedron ( id, i, j + 1, k ).y;
points[8] = blockgrid->Give_coord_hexahedron ( id, i, j + 1, k ).z;
points[20] = blockgrid->Give_coord_hexahedron ( id, i, j + 1, k + 1 ).z;
points[9] = blockgrid->Give_coord_hexahedron ( id, i + 1, j + 1, k ).x;
points[10] = blockgrid->Give_coord_hexahedron ( id, i + 1, j + 1, k ).y;
points[11] = blockgrid->Give_coord_hexahedron ( id, i + 1, j + 1, k ).z;
points[21] = blockgrid->Give_coord_hexahedron ( id, i + 1, j + 1, k + 1 ).x;
points[12] = blockgrid->Give_coord_hexahedron ( id, i, j, k + 1 ).x;
points[13] = blockgrid->Give_coord_hexahedron ( id, i, j, k + 1 ).y;
points[14] = blockgrid->Give_coord_hexahedron ( id, i, j, k + 1 ).z;
points[22] = blockgrid->Give_coord_hexahedron ( id, i + 1, j + 1, k + 1 ).y;
points[15] = blockgrid->Give_coord_hexahedron ( id, i + 1, j, k + 1 ).x;
points[16] = blockgrid->Give_coord_hexahedron ( id, i + 1, j, k + 1 ).y;
points[17] = blockgrid->Give_coord_hexahedron ( id, i + 1, j, k + 1 ).z;
points[23] = blockgrid->Give_coord_hexahedron ( id, i + 1, j + 1, k + 1 ).z;
points[18] = blockgrid->Give_coord_hexahedron ( id, i, j + 1, k + 1 ).x;
points[19] = blockgrid->Give_coord_hexahedron ( id, i, j + 1, k + 1 ).y;
points[20] = blockgrid->Give_coord_hexahedron ( id, i, j + 1, k + 1 ).z;
points[21] = blockgrid->Give_coord_hexahedron ( id, i + 1, j + 1, k + 1 ).x;
points[22] = blockgrid->Give_coord_hexahedron ( id, i + 1, j + 1, k + 1 ).y;
points[23] = blockgrid->Give_coord_hexahedron ( id, i + 1, j + 1, k + 1 ).z;
G ( points );
G(points);
loc_stiff = G.integrate ( bilinear_form );
loc_stiff = G.integrate ( bilinear_form );
for(int v = 0;v < 8;++v )
for(int u = 0;u < 8;++u )
loc_m_hexahedra[id][Ind_loc_matrix_hexahedra ( i,j,k ) *64+8*v+u]
for(int v = 0;v < 8;++v )
for(int u = 0;u < 8;++u )
loc_m_hexahedra[id][Ind_loc_matrix_hexahedra ( i,j,k ) *64+8*v+u]
= loc_stiff[u][v];
}
}
}
}
}
}
}
}
Update();
}
......
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