Neue funktion Calculate_boundary hinzugefügt, welche einen Marker übernimmt....

Neue funktion Calculate_boundary hinzugefügt, welche einen Marker übernimmt. so kann angegeben werden, für welchen rand die LSM berechnet werden soll.

program/source/grid/blockgrid.cc

@@ -87,9 +87,15 @@ void Blockgrid::init_blockgrid_coordinates()
         int Ny = this->Give_Ny_hexahedron(id_hex)+1;
         int Nz = this->Give_Nz_hexahedron(id_hex)+1;
 
-            for(int i=0;i<Nx;++i) for(int j=0;j<Nz;++j) for(int k=0;k<Nz;++k)
+            for(int i=0;i<Nx;++i) for(int j=0;j<Ny;++j) for(int k=0;k<Nz;++k)
             {
+                D3vector test(0,0,0);
+
                 blockgrid_hexa_coordinates.at(id_hex).at(    i    +Nx*(j    +Ny* k)) = Give_coord_hexahedron(id_hex,i,j,k);
+                if (blockgrid_hexa_coordinates.at(id_hex).at(    i    +Nx*(j    +Ny* k)) == test)
+                {
+                    cout << "debug " << endl;
+                }
                 blockgrid_hexa_coordinates_set.at(id_hex).at(i    +Nx*(j    +Ny* k)) = true;
             }
     }

program/source/grid/blockgrid.h

@@ -77,7 +77,14 @@ class  Blockgrid {
 
   int Total_number_of_points() const;
 
+
+  /*
+   * berechnet werte an allen knoten und speichert diese. Beschleunigt Give_coord_**** Funktionen deutlich.
+   * */
   void init_blockgrid_coordinates();
+  /*
+   * berechnet indizies der randkoordinaten, sodass direkt die benachabrten konten am Rand von hexahedrons angesprochen werden können.
+   * */
   void init_blockgrid_coordinates_boundary();
   void test_blockgrid_coordinates();
   void delete_blockgrid_coordinates();

program/source/grid/examples_ug_optics.cc

@@ -872,9 +872,12 @@ Lens::Lens (double Radius, double radius, double thickness, double curvatureLeft
   construction_done();
 }
 
-LensQuadrangle::LensQuadrangle(double Radius, double radius, double thickness, double curvatureLeft, double curvatureRight, double offsetX, double offsetY, double offsetZ, bool inner_grid_arched)
+LensQuadrangle::LensQuadrangle(double Radius, double thickness, double curvatureLeft, double curvatureRight, double offsetX, double offsetY, double offsetZ, bool inner_grid_arched, double radius)
 {
-
+    if (radius == 0)
+    {
+        radius = Radius / 2.0;
+    }
     assert(Radius > radius);
     size_pointer_global_data = 10;
     if (fabs(curvatureLeft) >1e10)

program/source/grid/examples_ug_optics.h

@@ -21,7 +21,7 @@ class Lens : public Unstructured_grid {
 
 class LensQuadrangle : public Unstructured_grid {
     public:
-    LensQuadrangle(double Radius, double radius, double thickness, double curvatureLeft, double curvatureRight, double offsetX, double offsetY, double offsetZ, bool inner_grid_arched = false);
+    LensQuadrangle(double Radius, double thickness, double curvatureLeft, double curvatureRight, double offsetX, double offsetY, double offsetZ, bool inner_grid_arched = false, double radius = 0.0);
     ~LensQuadrangle(){};
 
     double Radius;

program/source/grid/marker.cc

@@ -176,9 +176,10 @@ void Marker::selectHexBlocks ( unsigned int firstBlock, unsigned int numBlocks)
 	this->off();
 	for ( unsigned int i = 0; i < numBlocks; ++i ) {
 		selectHexBlock(firstBlock+i);
-	}
+    }
 }
 
+
 void Marker::selectHexBlocks ( std::vector<unsigned int> & chosenHexBlocks ) {
 	this->off();
 	for ( unsigned int count = 0;
@@ -189,6 +190,28 @@ void Marker::selectHexBlocks ( std::vector<unsigned int> & chosenHexBlocks ) {
 	}
 }
 
+void Marker::selectQuad(unsigned int chosenQuad)
+{
+    if ( ( int ) chosenQuad >= ug->Give_number_quadrangles() )
+        throw ( std::runtime_error ( std::string ( "Marker::selectQuad(): Selected quad does not exist" ) ) );
+    Set_marker<quadrangleEl> ( chosenQuad, yes_mark );
+
+
+}
+
+void Marker::selectSingleQuad(unsigned int chosenQuad)
+{
+    this->off();
+    selectQuad(chosenQuad);
+}
+
+void Marker::selectQuads(unsigned int firstQuad, unsigned int numQuads)
+{
+    this->off();
+    for ( unsigned int i = 0; i < numQuads; ++i ) {
+        selectQuad(firstQuad+i);
+    }
+}
 
 Marker::~Marker() {
 	if ( marker_hexahedra != NULL )
@@ -232,6 +255,25 @@ Boundary_Marker::Boundary_Marker ( Unstructured_grid* ug_ ) : Marker ( ug_ )
 
 Boundary_Marker::~Boundary_Marker() {}
 
+void Boundary_Marker::markBoundaryDir3D(dir3D dir)
+{
+    Quadrangle_el* quad;
+    Hexahedron_el* hex ;
+    for ( int id = 0;id < ug->Give_number_hexahedra();++id ) {
+            hex = ug->Give_hexahedron ( id );
+
+            if ( ug->Give_quadrangle ( hex->Give_id_quadrangle ( dir ) )->Give_exists_exterior() == false ) {
+
+            quad = ug->Give_quadrangle ( hex->Give_id_quadrangle ( dir ) );
+            Set_marker<quadrangleEl> ( hex->Give_id_quadrangle ( dir ), yes_mark );
+            for ( int i=0; i<4; ++i )
+                Set_marker<edgeEl> ( quad->Give_id_edge ( ( dir2D ) i ), yes_mark );
+            for ( int i=0; i<4; ++i )
+                Set_marker<pointEl> ( quad->Give_id_corner ( ( dir2D_sons ) i ), yes_mark );
+        }
+    }
+}
+
 
 Unstructured_grid_Marker::Unstructured_grid_Marker ( Unstructured_grid* ug_ ) : Marker ( ug_ ) {
 	for ( int id=0; id<ug->Give_number_hexahedra(); ++id ) {

program/source/grid/marker.h

@@ -26,18 +26,25 @@ public:
 	Marker( Unstructured_grid* ug_ );
 	virtual ~Marker();
 
-	// select everything
+    //* select everything
 	void on();
 
-	//select nothing
+    //* select nothing
 	void off();
 
-	// select block with id chosenBlock
+    //*  select block (inkl quad, edge and corner elements) with id chosenBlock
 	void selectHexBlock(unsigned int chosenHexBlock);
+    //*  select block (inkl quad, edge and corner elements) with id chosenBlock and deselect all others
 	void selectSingleHexBlock ( unsigned int chosenHexBlock );
+    //*  select range of blocks (inkl quad, edge and corner elements)
 	void selectHexBlocks ( std::vector<unsigned int> & chosenHexBlocks );
 	void selectHexBlocks ( unsigned int firstBlock, unsigned int numBlocks);
 
+    //*  select quadrangle (excl edge and corner) with id chosenQuad
+    void selectQuad(unsigned int chosenQuad);
+    void selectSingleQuad(unsigned int chosenQuad);
+    void selectQuads( unsigned int firstQuad, unsigned int numQuads);
+
 	template <elementTyp TYP_EL>
 	inline marker_typ Give_marker(int i) const;
 
@@ -78,6 +85,8 @@ class Boundary_Marker : public Marker {
 public:
 	Boundary_Marker(Unstructured_grid* ug_);
 	~Boundary_Marker();
+    //* mark boundary in dirction dir3D (Ndir3D, Edir3D... )
+    void markBoundaryDir3D(dir3D dir);
 
 };
 

program/source/interpol/interpol.cc

@@ -46,10 +46,13 @@
 #include "../grid/examples_ug.h"
 #include "../grid/blockgrid.h"
 #include "../grid/marker.h"
-#include "../extemp/parallel.h"
 #include "../extemp/extemp.h"
+#include "../extemp/parallel.h"
+#include "../extemp/variable.h"
 #include "../extemp/co_fu.h"
 #include "../extemp/functor.h"
+
+
 #include "interpol.h"
 #include <sstream>
 #include <iomanip>

program/source/pde_op/diffbo.h

@@ -261,6 +261,250 @@ void Local_stiffness_matrix<TYPE2>::Calculate_boundary ( const B& bilinear_form
 }
 
 
+template <typename TYPE2>
+template <class B>
+void Local_stiffness_matrix<TYPE2>::Calculate_boundary ( const B& bilinear_form, Marker & marker, bool addToLocalStiffnesMatrix ) {
+    int ux, uy, uz, vx, vy, vz,   u, v;
+    int Nx, Ny, Nz, N_total;
+    Unstructured_grid *ug;
+    Hexahedron_el *hex;
+
+    std::vector<std::vector<TYPE2> > loc_stiff;
+    std::vector<double> points ( 12, 0. ); //Länge 12
+    _COLSAMM_::ELEMENTS::_Quadrangle3D_2<_COLSAMM_::Gauss2, TYPE2> G;
+    //_COLSAMM_::ELEMENTS::_Quadrangle3D_<_COLSAMM_::Gauss2,TYPE2> G;
+
+    ug = blockgrid->Give_unstructured_grid();
+
+    for ( int id = 0;id < ug->Give_number_hexahedra();++id ) {
+
+        if ( ug->Give_hexahedron ( id )->my_object ( my_rank ) ) {
+
+            hex = ug->Give_hexahedron ( id );
+
+            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 (!addToLocalStiffnesMatrix)
+            {
+                for ( int i = 0;i < Nx;++i ) for ( int j = 0;j < Ny;++j ) for ( int k = 0;k < Nz;++k ) {
+                            for ( v = 0;v < 8;++v ) for ( u = 0;u < 8;++u )
+                                    loc_m_hexahedra[id][Ind_loc_matrix_hexahedra ( i,j,k ) *64+8*v+u] = 0.0;
+                        }
+            }
+
+
+
+            if ( marker.Give_marker<quadrangleEl>( hex->Give_id_quadrangle ( (Wdir3D)) ) == yes_mark ) {
+                for ( int k = 0;k < Nz;++k )
+                    for ( int j = 0;j < Ny;++j ) {
+                        int i = 0;
+
+                        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, j + 1, k ).x;
+                        points[4] = blockgrid->Give_coord_hexahedron ( id, i, j + 1, k ).y;
+                        points[5] = blockgrid->Give_coord_hexahedron ( id, i, j + 1, k ).z;
+
+                        points[6] = blockgrid->Give_coord_hexahedron ( id, i, j, k + 1 ).x;
+                        points[7] = blockgrid->Give_coord_hexahedron ( id, i, j, k + 1 ).y;
+                        points[8] = blockgrid->Give_coord_hexahedron ( id, i, j, k + 1 ).z;
+
+                        points[9]  = blockgrid->Give_coord_hexahedron ( id, i, j + 1, k + 1 ).x;
+                        points[10] = blockgrid->Give_coord_hexahedron ( id, i, j + 1, k + 1 ).y;
+                        points[11] = blockgrid->Give_coord_hexahedron ( id, i, j + 1, k + 1 ).z;
+
+                        G ( points );
+                        loc_stiff = G.integrate ( bilinear_form );
+
+                        ux = vx = 0;
+
+                        for ( vy = 0;vy < 2;++vy ) for ( vz = 0;vz < 2;++vz )
+                                for ( uy = 0;uy < 2;++uy ) for ( uz = 0;uz < 2;++uz )
+                                        loc_m_hexahedra[id][Ind_loc_matrix_hexahedra ( i,j,k ) *64+8* ( 4*vz+2*vy+vx ) + ( 4*uz+2*uy+ux ) ]
+                                        += loc_stiff[2*uz+uy][2*vz+vy];
+                    }
+            }
+
+            if ( marker.Give_marker<quadrangleEl>( hex->Give_id_quadrangle ( (Edir3D)) ) == yes_mark ) {
+                for ( int k = 0;k < Nz;++k )
+                    for ( int j = 0;j < Ny;++j ) {
+                        int i = Nx - 1;
+
+                        points[0] = blockgrid->Give_coord_hexahedron ( id, i + 1, j, k ).x;
+                        points[1] = blockgrid->Give_coord_hexahedron ( id, i + 1, j, k ).y;
+                        points[2] = blockgrid->Give_coord_hexahedron ( id, i + 1, j, k ).z;
+
+                        points[3] = blockgrid->Give_coord_hexahedron ( id, i + 1, j + 1, k ).x;
+                        points[4] = blockgrid->Give_coord_hexahedron ( id, i + 1, j + 1, k ).y;
+                        points[5] = blockgrid->Give_coord_hexahedron ( id, i + 1, j + 1, k ).z;
+
+                        points[6] = blockgrid->Give_coord_hexahedron ( id, i + 1, j, k + 1 ).x;
+                        points[7] = blockgrid->Give_coord_hexahedron ( id, i + 1, j, k + 1 ).y;
+                        points[8] = blockgrid->Give_coord_hexahedron ( id, i + 1, j, k + 1 ).z;
+
+                        points[9]  = blockgrid->Give_coord_hexahedron ( id, i + 1, j + 1, k + 1 ).x;
+                        points[10] = blockgrid->Give_coord_hexahedron ( id, i + 1, j + 1, k + 1 ).y;
+                        points[11] = blockgrid->Give_coord_hexahedron ( id, i + 1, j + 1, k + 1 ).z;
+
+                        G ( points );
+                        loc_stiff = G.integrate ( bilinear_form );
+
+                        ux = vx = 1;
+
+                        for ( vy = 0;vy < 2;++vy ) for ( vz = 0;vz < 2;++vz )
+                                for ( uy = 0;uy < 2;++uy ) for ( uz = 0;uz < 2;++uz )
+                                        loc_m_hexahedra[id][Ind_loc_matrix_hexahedra ( i,j,k ) *64+8* ( 4*vz+2*vy+vx ) + ( 4*uz+2*uy+ux ) ]
+                                        += loc_stiff[2*uz+uy][2*vz+vy];
+                    }
+            }
+
+            if ( marker.Give_marker<quadrangleEl>( hex->Give_id_quadrangle ( (Sdir3D)) ) == yes_mark ) {
+                for ( int k = 0;k < Nz;++k )
+                    for ( int i = 0;i < Nx;++i ) {
+                        int j = 0;
+
+                        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, k + 1 ).x;
+                        points[7] = blockgrid->Give_coord_hexahedron ( id, i, j, k + 1 ).y;
+                        points[8] = blockgrid->Give_coord_hexahedron ( id, i, j, k + 1 ).z;
+
+                        points[9]  = blockgrid->Give_coord_hexahedron ( id, i + 1, j, k + 1 ).x;
+                        points[10] = blockgrid->Give_coord_hexahedron ( id, i + 1, j, k + 1 ).y;
+                        points[11] = blockgrid->Give_coord_hexahedron ( id, i + 1, j, k + 1 ).z;
+
+                        G ( points );
+                        loc_stiff = G.integrate ( bilinear_form );
+
+                        uy = vy = 0;
+
+                        for ( vx = 0;vx < 2;++vx ) for ( vz = 0;vz < 2;++vz )
+                                for ( ux = 0;ux < 2;++ux ) for ( uz = 0;uz < 2;++uz )
+                                        loc_m_hexahedra[id][Ind_loc_matrix_hexahedra ( i,j,k ) *64+8* ( 4*vz+2*vy+vx ) + ( 4*uz+2*uy+ux ) ]
+                                        += loc_stiff[2*uz+ux][2*vz+vx];
+                    }
+            }
+
+            if ( marker.Give_marker<quadrangleEl>( hex->Give_id_quadrangle ( (Ndir3D)) ) == yes_mark ) {
+                for ( int k = 0;k < Nz;++k )
+                    for ( int i = 0;i < Nx;++i ) {
+                        int j = Ny - 1;
+
+                        points[0] = blockgrid->Give_coord_hexahedron ( id, i, j + 1, k ).x;
+                        points[1] = blockgrid->Give_coord_hexahedron ( id, i, j + 1, k ).y;
+                        points[2] = blockgrid->Give_coord_hexahedron ( id, i, j + 1, k ).z;
+
+                        points[3] = blockgrid->Give_coord_hexahedron ( id, i + 1, j + 1, k ).x;
+                        points[4] = blockgrid->Give_coord_hexahedron ( id, i + 1, j + 1, k ).y;
+                        points[5] = blockgrid->Give_coord_hexahedron ( id, i + 1, j + 1, k ).z;
+
+                        points[6] = blockgrid->Give_coord_hexahedron ( id, i, j + 1, k + 1 ).x;
+                        points[7] = blockgrid->Give_coord_hexahedron ( id, i, j + 1, k + 1 ).y;
+                        points[8] = blockgrid->Give_coord_hexahedron ( id, i, j + 1, k + 1 ).z;
+
+                        points[9]  = blockgrid->Give_coord_hexahedron ( id, i + 1, j + 1, k + 1 ).x;
+                        points[10] = blockgrid->Give_coord_hexahedron ( id, i + 1, j + 1, k + 1 ).y;
+                        points[11] = blockgrid->Give_coord_hexahedron ( id, i + 1, j + 1, k + 1 ).z;
+
+                        G ( points );
+                        loc_stiff = G.integrate ( bilinear_form );
+
+                        uy = vy = 1;
+
+                        for ( vx = 0;vx < 2;++vx ) for ( vz = 0;vz < 2;++vz )
+                                for ( ux = 0;ux < 2;++ux ) for ( uz = 0;uz < 2;++uz )
+                                        loc_m_hexahedra[id][Ind_loc_matrix_hexahedra ( i,j,k ) *64+8* ( 4*vz+2*vy+vx ) + ( 4*uz+2*uy+ux ) ]
+                                        += loc_stiff[2*uz+ux][2*vz+vx];
+                    }
+            }
+
+            if ( marker.Give_marker<quadrangleEl>( hex->Give_id_quadrangle ( (Ddir3D)) ) == yes_mark ) {
+                for ( int j = 0;j < Ny;++j )
+                    for ( int i = 0;i < Nx;++i ) {
+                        int k = 0;
+
+                        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[8] = blockgrid->Give_coord_hexahedron ( id, i, j + 1, k ).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;
+
+                        G ( points );
+                        loc_stiff = G.integrate ( bilinear_form );
+
+                        uz = vz = 0;
+
+                        for ( vx = 0;vx < 2;++vx ) for ( vy = 0;vy < 2;++vy )
+                                for ( ux = 0;ux < 2;++ux ) for ( uy = 0;uy < 2;++uy )
+                                        loc_m_hexahedra[id][Ind_loc_matrix_hexahedra ( i,j,k ) *64+8* ( 4*vz+2*vy+vx ) + ( 4*uz+2*uy+ux ) ]
+                                        += loc_stiff[2*uy+ux][2*vy+vx];
+                    }
+            }
+
+            if ( marker.Give_marker<quadrangleEl>( hex->Give_id_quadrangle ( (Tdir3D)) ) == yes_mark ) {
+                for ( int j = 0;j < Ny;++j )
+                    for ( int i = 0;i < Nx;++i ) {
+                        int k = Nz - 1;
+
+                        points[0] = blockgrid->Give_coord_hexahedron ( id, i, j, k + 1 ).x;
+                        points[1] = blockgrid->Give_coord_hexahedron ( id, i, j, k + 1 ).y;
+                        points[2] = blockgrid->Give_coord_hexahedron ( id, i, j, k + 1 ).z;
+
+                        points[3] = blockgrid->Give_coord_hexahedron ( id, i + 1, j, k + 1 ).x;
+                        points[4] = blockgrid->Give_coord_hexahedron ( id, i + 1, j, k + 1 ).y;
+                        points[5] = blockgrid->Give_coord_hexahedron ( id, i + 1, j, k + 1 ).z;
+
+                        points[6] = blockgrid->Give_coord_hexahedron ( id, i, j + 1, k + 1 ).x;
+                        points[7] = blockgrid->Give_coord_hexahedron ( id, i, j + 1, k + 1 ).y;
+                        points[8] = blockgrid->Give_coord_hexahedron ( id, i, j + 1, k + 1 ).z;
+
+                        points[9]  = blockgrid->Give_coord_hexahedron ( id, i + 1, j + 1, k + 1 ).x;
+                        points[10] = blockgrid->Give_coord_hexahedron ( id, i + 1, j + 1, k + 1 ).y;
+                        points[11] = blockgrid->Give_coord_hexahedron ( id, i + 1, j + 1, k + 1 ).z;
+
+                        G ( points );
+                        loc_stiff = G.integrate ( bilinear_form );
+
+                        uz = vz = 1;
+
+                        for ( vx = 0;vx < 2;++vx ) for ( vy = 0;vy < 2;++vy )
+                                for ( ux = 0;ux < 2;++ux ) for ( uy = 0;uy < 2;++uy )
+                                        loc_m_hexahedra[id][Ind_loc_matrix_hexahedra ( i,j,k ) *64+8* ( 4*vz+2*vy+vx ) + ( 4*uz+2*uy+ux ) ]
+                                        += loc_stiff[2*uy+ux][2*vy+vx];
+                    }
+            }
+
+
+
+
+        }
+    }
+
+    Update();
+}
+
+
 // Zum Testen und Spielen von Matthias eingefuegt, integriert alle rechten (in positiver z-richtung (Top)) Endflächen
 
 template <typename TYPE2>

program/source/pde_op/diffop.h

@@ -644,11 +644,17 @@ class Local_stiffness_matrix : public Expr<Local_stiffness_matrix<TYPE2> > {
 		template <class B>
 		void Calculate_boundary ( const B& bilinear_form );
 		
+        template <class B>
+        void Calculate_boundary ( const B& bilinear_form, Marker & marker, bool addToLocalStiffnesMatrix );
+
 		template <class B>
 		void Calculate_all_right_boundaries ( const B& bilinear_form );
 		
 		template <class B>
 		void Calculate_all_left_boundaries(const B& bilinear_form);
+
+        template <class B>
+        void Calculate_all_radial_boundaries(const B& bilinear_form);
 		
 		
 		void Print();