From 608218d1fee686d5d0611d8b9c8d6470f5ddf4bb Mon Sep 17 00:00:00 2001 From: Wenqing Wang Date: Fri, 8 Jun 2018 14:09:05 +0200 Subject: [PATCH 1/2] Re-organized some includes --- FEM/fem_ele.cpp | 4 ++++ FEM/fem_ele.h | 11 +++-------- FEM/fem_ele_std.cpp | 4 ++++ FEM/fem_ele_vec.cpp | 5 +++++ FEM/mathlib.cpp | 7 +++++++ FEM/matrix_class.cpp | 3 ++- FEM/par_ddc.cpp | 25 +++++++++++++++++++------ FEM/par_ddc.h | 37 ++++++++++++++++++++----------------- FEM/pcs_dm.cpp | 4 ++++ FEM/rf_pcs.cpp | 4 ++++ 10 files changed, 72 insertions(+), 32 deletions(-) diff --git a/FEM/fem_ele.cpp b/FEM/fem_ele.cpp index 2a1aa56ad..9bd0edab8 100644 --- a/FEM/fem_ele.cpp +++ b/FEM/fem_ele.cpp @@ -22,6 +22,10 @@ #include "femlib.h" #include "mathlib.h" +#ifndef USE_PETSC +#include "par_ddc.h" +#endif + #include "ShapeFunctionPool.h" namespace FiniteElement diff --git a/FEM/fem_ele.h b/FEM/fem_ele.h index 2adee1262..4db445301 100644 --- a/FEM/fem_ele.h +++ b/FEM/fem_ele.h @@ -17,15 +17,10 @@ #define fem_INC -// C++ -//#include -//#include - -#if defined(USE_PETSC) // || defined(other parallel libs)//03.3012. WW #include "prototyp.h" -#else -// MSH -#include "par_ddc.h" //OK //Moved from fem_ele_std.h. WW + +#ifndef USE_PETSC +class CPARDomain; #endif #include "MSHEnums.h" diff --git a/FEM/fem_ele_std.cpp b/FEM/fem_ele_std.cpp index 324159179..198d8b816 100644 --- a/FEM/fem_ele_std.cpp +++ b/FEM/fem_ele_std.cpp @@ -50,6 +50,10 @@ using Math_Group::CSparseMatrix; #endif +#ifndef USE_PETSC +#include "par_ddc.h" +#endif + #ifdef OGS_USE_CVODE extern "C" { #include /* prototypes for CVODE fcts., consts. */ diff --git a/FEM/fem_ele_vec.cpp b/FEM/fem_ele_vec.cpp index 7f34fc9fc..5e5d3a7db 100644 --- a/FEM/fem_ele_vec.cpp +++ b/FEM/fem_ele_vec.cpp @@ -35,6 +35,11 @@ #include "equation_class.h" using Math_Group::CSparseMatrix; #endif + +#ifndef USE_PETSC +#include "par_ddc.h" +#endif + // #define COMP_MOL_MASS_AIR 28.96 // kg/kmol WW 28.96 #define GAS_CONSTANT 8314.41 // J/(kmol*K) WW diff --git a/FEM/mathlib.cpp b/FEM/mathlib.cpp index 1675759f2..ed914591a 100644 --- a/FEM/mathlib.cpp +++ b/FEM/mathlib.cpp @@ -161,15 +161,22 @@ /* Preprozessor-Definitionen */ #include +#include + //#include #include "femlib.h" //CMCD 03 2004 #include "makros.h" #include "memory.h" #include "display.h" #include "mathlib.h" +#include "geo_mathlib.h" // for M3KreuzProdukt + // WW---------------------- #include "par_ddc.h" // WW---------------------- + +#include "prototyp.h" + double pai = 4.0 * atan(1.0); VoidFuncDXCDX ShapeFunction; VoidFuncDXCDX ShapeFunctionHQ; diff --git a/FEM/matrix_class.cpp b/FEM/matrix_class.cpp index bd0d7751c..f0e6639b2 100644 --- a/FEM/matrix_class.cpp +++ b/FEM/matrix_class.cpp @@ -679,7 +679,8 @@ void vec::operator=(const vec& v) 03/2010 WW: CRS storage for matrix algbraic ******************************************************************** */ -SparseTable::SparseTable(CFEMesh* a_mesh, bool quadratic, bool symm, StorageType stype) +SparseTable::SparseTable(MeshLib::CFEMesh* a_mesh, bool quadratic, + bool symm, StorageType stype) : symmetry(symm), storage_type(stype) { long i = 0, j = 0, ii = 0, jj = 0; diff --git a/FEM/par_ddc.cpp b/FEM/par_ddc.cpp index 5b5565791..8a3e846cf 100644 --- a/FEM/par_ddc.cpp +++ b/FEM/par_ddc.cpp @@ -26,23 +26,24 @@ // #error to indicate a fatal error. Users can either #undef // the names before including mpi.h or include mpi.h *before* stdio.h // or iostream. +#include "par_ddc.h" + #if defined(USE_MPI) || defined(USE_MPI_PARPROC) || defined(USE_MPI_REGSOIL) || defined(USE_MPI_GEMS) \ || defined(USE_MPI_KRC) -//#undef SEEK_SET //WW -//#undef SEEK_END //WW -//#undef SEEK_CUR //WW #include #include "SplitMPI_Communicator.h" char t_fname[3]; double time_ele_paral; #endif -//---- MPI Parallel -------------- #include // C++ STL #include -using namespace std; -#include "par_ddc.h" + +#include "rf_pcs.h" + +#include "msh_mesh.h" + // FEM-Makros #include "makros.h" #ifndef NEW_EQS // WW. 11.2008 @@ -58,6 +59,18 @@ using namespace std; vector dom_vector; vector node_connected_doms; // This will be removed after sparse class is finished WW +using namespace std; +using process::CRFProcessDeformation; +using FiniteElement::CFiniteElementVec; + +#ifdef NEW_EQS +using Math_Group::Linear_EQS; +using Math_Group::SparseTable; +#endif + +#if defined(USE_MPI) +using MeshLib::CFEMesh; +#endif /************************************************************************** STRLib-Method: Task: diff --git a/FEM/par_ddc.h b/FEM/par_ddc.h index 1264898d0..ee553e394 100644 --- a/FEM/par_ddc.h +++ b/FEM/par_ddc.h @@ -24,7 +24,9 @@ #include #include -#include "rf_pcs.h" +#if defined(USE_MPI) +#include +#endif //--------- Declaration ---------------------WW namespace FiniteElement @@ -32,30 +34,31 @@ namespace FiniteElement class CFiniteElementStd; class CFiniteElementVec; } + +class CRFProcess; namespace process { class CRFProcessDeformation; } + +#ifdef NEW_EQS namespace Math_Group { class SparseTable; class Linear_EQS; } -using process::CRFProcessDeformation; -using FiniteElement::CFiniteElementVec; -using Math_Group::Linear_EQS; -using Math_Group::SparseTable; -//--------- Declaration ---------------------WW +class CNumerics; +#endif + +struct LINEAR_SOLVER; +struct LINEAR_SOLVER_PROPERTIES; -#if defined(USE_MPI) -// WW namespace MeshLib { class CFEMesh; } -using MeshLib::CFEMesh; -#endif -void FindNodesOnInterface(CFEMesh* m_msh, bool quadr); + +void FindNodesOnInterface(MeshLib::CFEMesh* m_msh, bool quadr); //----------------------------------------------- class CPARDomain @@ -69,7 +72,7 @@ class CPARDomain long num_inner_nodesHQ; long num_boundary_nodes; long num_boundary_nodesHQ; - friend void FindNodesOnInterface(CFEMesh* m_msh, bool quadr); + friend void FindNodesOnInterface(MeshLib::CFEMesh* m_msh, bool quadr); //#endif #if defined(USE_MPI) // 13.12.2007 WW // Store global indices of all border nodes to border_nodes of the whole mesh @@ -103,10 +106,10 @@ class CPARDomain long shift[5]; // WW // Equation #ifdef NEW_EQS // WW - SparseTable* sparse_graph; - SparseTable* sparse_graph_H; - Linear_EQS* eqs; // WW - Linear_EQS* eqsH; // WW + Math_Group::SparseTable* sparse_graph; + Math_Group::SparseTable* sparse_graph_H; + Math_Group::Linear_EQS* eqs; // WW + Math_Group::Linear_EQS* eqsH; // WW #endif friend class CRFProcess; // WW // WW //:: for SXC compiler @@ -128,7 +131,7 @@ class CPARDomain char* lsp_name; #endif // MSH - CFEMesh* m_msh; + MeshLib::CFEMesh* m_msh; // public: CPARDomain(void); ~CPARDomain(void); diff --git a/FEM/pcs_dm.cpp b/FEM/pcs_dm.cpp index 53399220f..3c4539f6c 100644 --- a/FEM/pcs_dm.cpp +++ b/FEM/pcs_dm.cpp @@ -54,6 +54,10 @@ #include "PETSC/PETScLinearSolver.h" #endif +#ifndef USE_PETSC +#include "par_ddc.h" +#endif + using namespace std; // Solver diff --git a/FEM/rf_pcs.cpp b/FEM/rf_pcs.cpp index 47753139b..d0613597c 100644 --- a/FEM/rf_pcs.cpp +++ b/FEM/rf_pcs.cpp @@ -130,6 +130,10 @@ REACT_BRNS* m_vec_BRNS; #include "fct_mpi.h" +#ifndef USE_PETSC +#include "par_ddc.h" +#endif + using namespace std; using namespace MeshLib; using namespace Math_Group; From 211ee226780a43a344594356058ca4918bb3e0c5 Mon Sep 17 00:00:00 2001 From: Wenqing Wang Date: Fri, 8 Jun 2018 14:12:50 +0200 Subject: [PATCH 2/2] Removed pieces of unused code --- FEM/mathlib.cpp | 2498 ----------------------------------------------- FEM/mathlib.h | 124 --- 2 files changed, 2622 deletions(-) diff --git a/FEM/mathlib.cpp b/FEM/mathlib.cpp index ed914591a..5a1533204 100644 --- a/FEM/mathlib.cpp +++ b/FEM/mathlib.cpp @@ -210,2504 +210,6 @@ double MBtrgVec(double* vec, long n) return sqrt(zwo); } -#ifdef obsolete // 05.03.2010 WW -/*************************************************************************** - ROCKFLOW - Funktion: MGleichDouble - Aufgabe: - Vergleicht zwei double-Zahlen unter Beruecksichtigung - einer Fehlertoleranz - Formalparameter: - E: zahl1, zahl2 - E: tol - Fehlertoleranz (positiv) - Ergebnis: - 0 - ungleich - 1 - gleich - Programmaenderungen: - 07/1994 hh Erste Version - - **************************************************************************/ - -int MGleichDouble(double zahl1, double zahl2, double tol) -{ - int retval; - if (fabs(zahl1 - zahl2) <= tol) - retval = 1; - else - retval = 0; - return retval; -} -#endif //#ifdef obsolete //05.03.2010 WW - -//////////////////////////////////////////////////////////// -#ifdef obsolete // 05.03.2010 WW -/*************************************************************************** - ROCKFLOW - Funktion: MOmega1D - Aufgabe: - Berechnet Ansatzfunktion (r). - / \ - 1 | (1-r) | - Omega = --- | | - 2 | (1+r) | - \ / - Formalparameter: - Z: *vf - 1x2 Feld - E: r (s.o.) - Ergebnis: - Vektor - Aenderungen/Korrekturen: - 08/1995 cb Erste Version - - **************************************************************************/ - -int MOmega1D(double* vf, double r) -{ - long i; - int ok = 0; -#ifdef ERROR_CONTROL - if (vf == NULL) - { - ok = 0; - return ok; - } -#endif - vf[0] = (1.0 + r); - vf[1] = (1.0 - r); - for (i = 0; i < 2; i++) - vf[i] *= 0.5; - return ok = 1; -} - -/*************************************************************************** - ROCKFLOW - Funktion: MOmega2D - Aufgabe: - Berechnet Ansatzfunktion (r,s). - / (1+r)(1+s) \ - | | - 1 | (1-r)(1+s) | - Omega = --- | | - 4 | (1-r)(1-s) | - | | - \ (1+r)(1-s) / - Formalparameter: - Z: *vf - 1x4 Feld - E: r,s (s.o.) - Ergebnis: - Vektor - Aenderungen/Korrekturen: - 08/1995 cb Erste Version - - **************************************************************************/ - -int MOmega2D(double* vf, double r, double s) -{ - long i; - int ok = 0; -#ifdef ERROR_CONTROL - if (vf == NULL) - { - ok = 0; - return ok; - } -#endif - vf[0] = (1.0 + r) * (1.0 + s); - vf[1] = (1.0 - r) * (1.0 + s); - vf[2] = (1.0 - r) * (1.0 - s); - vf[3] = (1.0 + r) * (1.0 - s); - for (i = 0; i < 4; i++) - vf[i] *= 0.25; - return ok = 1; -} - -/*************************************************************************** - ROCKFLOW - Funktion: MOmega3D - Aufgabe: - Berechnet Ansatzfunktion (r,s). - / (1+r)(1+s)(1+t) \ - | (1-r)(1+s)(1+t) | - 1 | (1-r)(1-s)(1+t) | - Omega = --- | (1+r)(1-s)(1+t) | - 8 | (1+r)(1+s)(1-t) | - | (1-r)(1+s)(1-t) | - | (1-r)(1-s)(1-t) | - \ (1+r)(1-s)(1-t) / - Formalparameter: - Z: *vf - 1x8 Feld - E: r,s,t (s.o.) - Ergebnis: - Vektor - Aenderungen/Korrekturen: - 08/1995 cb Erste Version - - **************************************************************************/ - -int MOmega3D(double* vf, double r, double s, double t) -{ - long i; - int ok = 0; -#ifdef ERROR_CONTROL - if (vf == NULL) - { - ok = 0; - return ok; - } -#endif - vf[0] = (1.0 + r) * (1.0 + s) * (1.0 + t); - vf[1] = (1.0 - r) * (1.0 + s) * (1.0 + t); - vf[2] = (1.0 - r) * (1.0 - s) * (1.0 + t); - vf[3] = (1.0 + r) * (1.0 - s) * (1.0 + t); - vf[4] = (1.0 + r) * (1.0 + s) * (1.0 - t); - vf[5] = (1.0 - r) * (1.0 + s) * (1.0 - t); - vf[6] = (1.0 - r) * (1.0 - s) * (1.0 - t); - vf[7] = (1.0 + r) * (1.0 - s) * (1.0 - t); - for (i = 0; i < 8; i++) - vf[i] *= 0.125; - return ok = 1; -} - -/*************************************************************************** - ROCKFLOW - Funktion: MOmega2DTriangle - Aufgabe: - Berechnet Ansatzfunktion (r,s). - / (1+r)(1+s) \ - | | - 1 | (1-r)(1+s) | - Omega = --- | | - 4 | (1-r)(1-s) | - | | - \ (1+r)(1-s) / - Formalparameter: - Z: *vf - 1x4 Feld - E: r,s (s.o.) - Ergebnis: - Vektor - Aenderungen/Korrekturen: - 08/1995 cb Erste Version - - **************************************************************************/ -//#ifdef obsolete //WW. 06.11.2008 -int MOmega2DTriangle(double* vf, double xx, double yy, long number) -{ - int i, nn = 3, ok = 0; - double x[3], y[3]; - double area; - long* element_nodes; - -#ifdef ERROR_CONTROL - if (vf == NULL) - { - ok = 0; - return ok; - } -#endif - area = fabs(ElGetElementVolume(number)); // should be triangle area ! - /* Calc2DElementCoordinatesTriangle(number,x,y); */ - element_nodes = ElGetElementNodes(number); - for (i = 0; i < nn; i++) - { - x[i] = GetNodeX(element_nodes[i]); - y[i] = GetNodeY(element_nodes[i]); - } - /* CMCD March 2004 The area of the triangle needs to be calculated for the 2D projected triangle - not the three D as above - area = (x[1]*y[2]-x[2]*y[1]+x[0]*y[1]-x[1]*y[0]+x[2]*y[0]-x[0]*y[2])/2; // CMCD March 2004 - xx=0.; - yy=0.; - what is this ? - */ - vf[0] = (x[1] * y[2] - x[2] * y[1]) + (y[1] - y[2]) * xx + (x[2] - x[1]) * yy; - vf[1] = (x[2] * y[0] - x[0] * y[2]) + (y[2] - y[0]) * xx + (x[0] - x[2]) * yy; - vf[2] = (x[0] * y[1] - x[1] * y[0]) + (y[0] - y[1]) * xx + (x[1] - x[0]) * yy; - vf[0] /= (2. * area); - vf[1] /= (2. * area); - vf[2] /= (2. * area); - - return ok = 1; -} - -/************************************************************************** - GeoLib-Method: MOmega3DTetrahedron - Task: Interpolation function for tetrahedra - Programing: - 08/2003 OK Implementation -**************************************************************************/ -int MOmega3DTetrahedron(double* vf, double r, double s, double t, long number) -{ - int ok = 0; - long* element_nodes; - int i; - double a1, a2, a3, a4; - double b1, b2, b3, b4; - double c1, c2, c3, c4; - double d1, d2, d3, d4; - double x[4], y[4], z[4]; - double N1, N2, N3, N4; - double volume; - double mat3x3[9]; - - volume = fabs(ElGetElementVolume(number)); - element_nodes = ElGetElementNodes(number); - for (i = 0; i < 4; i++) - { - x[i] = GetNodeX(element_nodes[i]); - y[i] = GetNodeY(element_nodes[i]); - z[i] = GetNodeZ(element_nodes[i]); - } - - mat3x3[0] = x[1]; - mat3x3[1] = y[1]; - mat3x3[2] = z[1]; - mat3x3[3] = x[2]; - mat3x3[4] = y[2]; - mat3x3[5] = z[2]; - mat3x3[6] = x[3]; - mat3x3[7] = y[3]; - mat3x3[8] = z[3]; - a1 = -1.0 * M3Determinante(mat3x3); - mat3x3[0] = x[2]; - mat3x3[1] = y[2]; - mat3x3[2] = z[2]; - mat3x3[3] = x[3]; - mat3x3[4] = y[3]; - mat3x3[5] = z[3]; - mat3x3[6] = x[0]; - mat3x3[7] = y[0]; - mat3x3[8] = z[0]; - a2 = 1.0 * M3Determinante(mat3x3); - mat3x3[0] = x[3]; - mat3x3[1] = y[3]; - mat3x3[2] = z[3]; - mat3x3[3] = x[0]; - mat3x3[4] = y[0]; - mat3x3[5] = z[0]; - mat3x3[6] = x[1]; - mat3x3[7] = y[1]; - mat3x3[8] = z[1]; - a3 = -1.0 * M3Determinante(mat3x3); - mat3x3[0] = x[0]; - mat3x3[1] = y[0]; - mat3x3[2] = z[0]; - mat3x3[3] = x[1]; - mat3x3[4] = y[1]; - mat3x3[5] = z[1]; - mat3x3[6] = x[2]; - mat3x3[7] = y[2]; - mat3x3[8] = z[2]; - a4 = 1.0 * M3Determinante(mat3x3); - - mat3x3[0] = 1.0; - mat3x3[1] = y[1]; - mat3x3[2] = z[1]; - mat3x3[3] = 1.0; - mat3x3[4] = y[2]; - mat3x3[5] = z[2]; - mat3x3[6] = 1.0; - mat3x3[7] = y[3]; - mat3x3[8] = z[3]; - b1 = -1.0 * M3Determinante(mat3x3); - mat3x3[0] = 1.0; - mat3x3[1] = y[2]; - mat3x3[2] = z[2]; - mat3x3[3] = 1.0; - mat3x3[4] = y[3]; - mat3x3[5] = z[3]; - mat3x3[6] = 1.0; - mat3x3[7] = y[0]; - mat3x3[8] = z[0]; - b2 = 1.0 * M3Determinante(mat3x3); - mat3x3[0] = 1.0; - mat3x3[1] = y[3]; - mat3x3[2] = z[3]; - mat3x3[3] = 1.0; - mat3x3[4] = y[0]; - mat3x3[5] = z[0]; - mat3x3[6] = 1.0; - mat3x3[7] = y[1]; - mat3x3[8] = z[1]; - b3 = -1.0 * M3Determinante(mat3x3); - mat3x3[0] = 1.0; - mat3x3[1] = y[0]; - mat3x3[2] = z[0]; - mat3x3[3] = 1.0; - mat3x3[4] = y[1]; - mat3x3[5] = z[1]; - mat3x3[6] = 1.0; - mat3x3[7] = y[2]; - mat3x3[8] = z[2]; - b4 = 1.0 * M3Determinante(mat3x3); - - mat3x3[0] = x[1]; - mat3x3[1] = 1.0; - mat3x3[2] = z[1]; - mat3x3[3] = x[2]; - mat3x3[4] = 1.0; - mat3x3[5] = z[2]; - mat3x3[6] = x[3]; - mat3x3[7] = 1.0; - mat3x3[8] = z[3]; - c1 = -1.0 * M3Determinante(mat3x3); - mat3x3[0] = x[2]; - mat3x3[1] = 1.0; - mat3x3[2] = z[2]; - mat3x3[3] = x[3]; - mat3x3[4] = 1.0; - mat3x3[5] = z[3]; - mat3x3[6] = x[0]; - mat3x3[7] = 1.0; - mat3x3[8] = z[0]; - c2 = 1.0 * M3Determinante(mat3x3); - mat3x3[0] = x[3]; - mat3x3[1] = 1.0; - mat3x3[2] = z[3]; - mat3x3[3] = x[0]; - mat3x3[4] = 1.0; - mat3x3[5] = z[0]; - mat3x3[6] = x[1]; - mat3x3[7] = 1.0; - mat3x3[8] = z[1]; - c3 = -1.0 * M3Determinante(mat3x3); - mat3x3[0] = x[0]; - mat3x3[1] = 1.0; - mat3x3[2] = z[0]; - mat3x3[3] = x[1]; - mat3x3[4] = 1.0; - mat3x3[5] = z[1]; - mat3x3[6] = x[2]; - mat3x3[7] = 1.0; - mat3x3[8] = z[2]; - c4 = 1.0 * M3Determinante(mat3x3); - - mat3x3[0] = x[1]; - mat3x3[1] = y[1]; - mat3x3[2] = 1.0; - mat3x3[3] = x[2]; - mat3x3[4] = y[2]; - mat3x3[5] = 1.0; - mat3x3[6] = x[3]; - mat3x3[7] = y[3]; - mat3x3[8] = 1.0; - d1 = -1.0 * M3Determinante(mat3x3); - mat3x3[0] = x[2]; - mat3x3[1] = y[2]; - mat3x3[2] = 1.0; - mat3x3[3] = x[3]; - mat3x3[4] = y[3]; - mat3x3[5] = 1.0; - mat3x3[6] = x[0]; - mat3x3[7] = y[0]; - mat3x3[8] = 1.0; - d2 = 1.0 * M3Determinante(mat3x3); - mat3x3[0] = x[3]; - mat3x3[1] = y[3]; - mat3x3[2] = 1.0; - mat3x3[3] = x[0]; - mat3x3[4] = y[0]; - mat3x3[5] = 1.0; - mat3x3[6] = x[1]; - mat3x3[7] = y[1]; - mat3x3[8] = 1.0; - d3 = -1.0 * M3Determinante(mat3x3); - mat3x3[0] = x[0]; - mat3x3[1] = y[0]; - mat3x3[2] = 1.0; - mat3x3[3] = x[1]; - mat3x3[4] = y[1]; - mat3x3[5] = 1.0; - mat3x3[6] = x[2]; - mat3x3[7] = y[2]; - mat3x3[8] = 1.0; - d4 = 1.0 * M3Determinante(mat3x3); - - // Element Shape Functions - N1 = ((a1 * 1) + (b1 * r) + (c1 * s) + (d1 * t)) / (6 * volume); - N2 = ((a2 * 1) + (b2 * r) + (c2 * s) + (d2 * t)) / (6 * volume); - N3 = ((a3 * 1) + (b3 * r) + (c3 * s) + (d3 * t)) / (6 * volume); - N4 = ((a4 * 1) + (b4 * r) + (c4 * s) + (d4 * t)) / (6 * volume); - - vf[0] = N1; - vf[1] = N2; - vf[2] = N3; - vf[3] = N4; - - return ok = 1; -} - -//#endif //#ifndef obsolete //WW. 06.11.2008 - -/*************************************************************************** - ROCKFLOW - Funktion: MPhi3D - Aufgabe: - Berechnet Testfunktion (r,s,t) ohne Upwind-Parameter alpha. - Phi = Omega + alpha * grad Omega - Formalparameter: - Z: *vf - 1x8 Feld - E: r,s,t (s.o.) - Ergebnis: - Vektor - Aenderungen/Korrekturen: - 01.07.1997 R.Kaiser erste Version - - **************************************************************************/ - -int MPhi3D(double* vf, double r, double s, double t) -{ - int i; - int ok = 0; -#ifdef ERROR_CONTROL - if (vf == NULL) - { - ok = 0; - return ok; - } -#endif - vf[0] = (1.0 + r) * (1.0 + s) * (1.0 + t); - vf[1] = (1.0 - r) * (1.0 + s) * (1.0 + t); - vf[2] = (1.0 - r) * (1.0 - s) * (1.0 + t); - vf[3] = (1.0 + r) * (1.0 - s) * (1.0 + t); - vf[4] = (1.0 + r) * (1.0 + s) * (1.0 - t); - vf[5] = (1.0 - r) * (1.0 + s) * (1.0 - t); - vf[6] = (1.0 - r) * (1.0 - s) * (1.0 - t); - vf[7] = (1.0 + r) * (1.0 - s) * (1.0 - t); - for (i = 0; i < 8; i++) - vf[i] *= 0.125; - return ok = 1; -} - -/*************************************************************************** - ROCKFLOW - Funktion: MPhi2D_SUPG - Aufgabe: - Berechnet Testfunktion (r,s) incl. Upwind-Parameter alpha. - Phi = Omega + alpha * grad Omega - Formalparameter: - Z: *vf - 1x4 Feld - E: r,s (s.o.) - Ergebnis: - Vektor - Aenderungen/Korrekturen: - 11/1995 cb Erste Version - 01.07.1997 R.Kaiser Umbenannt zu MPhi2D_SUPG - - **************************************************************************/ - -int MPhi2D_SUPG(double* vf, double r, double s, double* alpha) -{ - int i; - int ok = 0; -#ifdef ERROR_CONTROL - if (vf == NULL) - { - ok = 0; - return ok; - } -#endif - vf[0] = (1.0 + r) * (1.0 + s) + alpha[0] * (1.0 + s) + alpha[1] * (1.0 + r); - vf[1] = (1.0 - r) * (1.0 + s) - alpha[0] * (1.0 + s) + alpha[1] * (1.0 - r); - vf[2] = (1.0 - r) * (1.0 - s) - alpha[0] * (1.0 - s) - alpha[1] * (1.0 - r); - vf[3] = (1.0 + r) * (1.0 - s) + alpha[0] * (1.0 - s) - alpha[1] * (1.0 + r); - for (i = 0; i < 4; i++) - vf[i] *= 0.25; - return ok = 1; -} - -/*************************************************************************** - ROCKFLOW - Funktion: MPhi3D_SUPG - Aufgabe: - Berechnet Testfunktion (r,s,t) incl. Upwind-Parameter alpha. - Phi = Omega + alpha * grad Omega - Formalparameter: - Z: *vf - 1x8 Feld - E: r,s,t (s.o.) - Ergebnis: - Vektor - Aenderungen/Korrekturen: - 11/1995 cb Erste Version - 01.07.1997 R.Kaiser Umbenannt zu MPhi3D_SUPG - - **************************************************************************/ - -int MPhi3D_SUPG(double* vf, double r, double s, double t, double* alpha) -{ - int i; - int ok = 0; -#ifdef ERROR_CONTROL - if (vf == NULL) - { - ok = 0; - return ok; - } -#endif - vf[0] = (1.0 + r) * (1.0 + s) * (1.0 + t) + alpha[0] * (1.0 + s) * (1.0 + t) + alpha[1] * (1.0 + r) * (1.0 + t) - + alpha[2] * (1.0 + r) * (1.0 + s); - vf[1] = (1.0 - r) * (1.0 + s) * (1.0 + t) - alpha[0] * (1.0 + s) * (1.0 + t) + alpha[1] * (1.0 - r) * (1.0 + t) - + alpha[2] * (1.0 - r) * (1.0 + s); - vf[2] = (1.0 - r) * (1.0 - s) * (1.0 + t) - alpha[0] * (1.0 - s) * (1.0 + t) - alpha[1] * (1.0 - r) * (1.0 + t) - + alpha[2] * (1.0 - r) * (1.0 - s); - vf[3] = (1.0 + r) * (1.0 - s) * (1.0 + t) + alpha[0] * (1.0 - s) * (1.0 + t) - alpha[1] * (1.0 + r) * (1.0 + t) - + alpha[2] * (1.0 + r) * (1.0 - s); - vf[4] = (1.0 + r) * (1.0 + s) * (1.0 - t) + alpha[0] * (1.0 + s) * (1.0 - t) + alpha[1] * (1.0 + r) * (1.0 - t) - - alpha[2] * (1.0 + r) * (1.0 + s); - vf[5] = (1.0 - r) * (1.0 + s) * (1.0 - t) - alpha[0] * (1.0 + s) * (1.0 - t) + alpha[1] * (1.0 - r) * (1.0 - t) - - alpha[2] * (1.0 - r) * (1.0 + s); - vf[6] = (1.0 - r) * (1.0 - s) * (1.0 - t) - alpha[0] * (1.0 - s) * (1.0 - t) - alpha[1] * (1.0 - r) * (1.0 - t) - - alpha[2] * (1.0 - r) * (1.0 - s); - vf[7] = (1.0 + r) * (1.0 - s) * (1.0 - t) + alpha[0] * (1.0 - s) * (1.0 - t) - alpha[1] * (1.0 + r) * (1.0 - t) - - alpha[2] * (1.0 + r) * (1.0 - s); - for (i = 0; i < 8; i++) - vf[i] *= 0.125; - return ok = 1; -} - -/*************************************************************************** - ROCKFLOW - Funktion: MGradOmega2D - Aufgabe: - Berechnet Gradient eines 2D Vektorfeldes, - dessen Ansatzfunktion (Vektor) bekannt ist (r,s). - / \ - 1 | +(1+s) -(1+s) -(1-s) +(1-s) | - grad (omega)= --- | | - 4 | +(1+r) +(1-r) -(1-r) -(1+r) | - \ / - Formalparameter: - Z: *vf - 2x4 Feld - E: r,s (s.o.) - Ergebnis: - 2x4 Matrix - Aenderungen/Korrekturen: - 05/1995 hh Erste Version - - **************************************************************************/ - -int MGradOmega2D(double* vf, double r, double s) -{ - long i; - int ok = 0; -#ifdef ERROR_CONTROL - if (vf == NULL) - { - ok = 0; - return ok; - } -#endif - vf[0] = +(1.0 + s); - vf[1] = -(1.0 + s); - vf[2] = -(1.0 - s); - vf[3] = +(1.0 - s); - vf[4] = +(1.0 + r); - vf[5] = +(1.0 - r); - vf[6] = -(1.0 - r); - vf[7] = -(1.0 + r); - for (i = 0; i < 8; i++) - vf[i] *= 0.25; - return ok = 1; -} - -/*************************************************************************** - ROCKFLOW - Funktion: MGradOmega3D - Aufgabe: - Berechnet Gradient eines 3D Vektorfeldes, - dessen Ansatzfunktion (Vektor) bekannt ist (r,s,t). - - | 0 1 2 3 4 5 6 7 | - | | - grad (omega)= | 8 9 10 11 12 13 14 15 | - | | - | 16 17 18 19 20 21 22 23 | - - Zahlen in der Matrix stehen fuer Positionen im Feld vf. - Gleichungen s.u.. - - Formalparameter: - Z: *vf - 3x8 Feld - E: r,s (s.o.) - Ergebnis: - 3x8 Matrix - Aenderungen/Korrekturen: - 05/1995 hh Erste Version - - **************************************************************************/ - -int MGradOmega3D(double* vf, double r, double s, double t) -{ - long i; - int ok = 0; -#ifdef ERROR_CONTROL - if (vf == NULL) - { - ok = 0; - return ok; - } -#endif - vf[0] = +(1.0 + s) * (1.0 + t); - vf[1] = -(1.0 + s) * (1.0 + t); - vf[2] = -(1.0 - s) * (1.0 + t); - vf[3] = +(1.0 - s) * (1.0 + t); - - vf[4] = +(1.0 + s) * (1.0 - t); - vf[5] = -(1.0 + s) * (1.0 - t); - vf[6] = -(1.0 - s) * (1.0 - t); - vf[7] = +(1.0 - s) * (1.0 - t); - - vf[8] = +(1.0 + r) * (1.0 + t); - vf[9] = +(1.0 - r) * (1.0 + t); - vf[10] = -(1.0 - r) * (1.0 + t); - vf[11] = -(1.0 + r) * (1.0 + t); - - vf[12] = +(1.0 + r) * (1.0 - t); - vf[13] = +(1.0 - r) * (1.0 - t); - vf[14] = -(1.0 - r) * (1.0 - t); - vf[15] = -(1.0 + r) * (1.0 - t); - - vf[16] = +(1.0 + r) * (1.0 + s); - vf[17] = +(1.0 - r) * (1.0 + s); - vf[18] = +(1.0 - r) * (1.0 - s); - vf[19] = +(1.0 + r) * (1.0 - s); - - vf[20] = -(1.0 + r) * (1.0 + s); - vf[21] = -(1.0 - r) * (1.0 + s); - vf[22] = -(1.0 - r) * (1.0 - s); - vf[23] = -(1.0 + r) * (1.0 - s); - - for (i = 0; i < 24; i++) - vf[i] *= 0.125; - return ok = 1; -} - -/*************************************************************************** - ROCKFLOW - Funktion: MGradPhi2D - Aufgabe: - Berechnet Gradient der Testfunktionen (r,s) - grad phi = grad omega - Formalparameter: - Z: *vf - 2x4 Feld - E: r,s (s.o.) - Ergebnis: - 2x4 Matrix - Aenderungen/Korrekturen: - 11/1995 cb Erste Version - - **************************************************************************/ - -int MGradPhi2D(double* vf, double r, double s) -{ - return MGradOmega2D(vf, r, s); -} - -/*************************************************************************** - ROCKFLOW - Funktion: MGradPhi3D - Aufgabe: - Berechnet Gradient der Testfunktionen (r,s,t). - grad phi = grad omega - Formalparameter: - Z: *vf - 3x8 Feld - E: r,s (s.o.) - Ergebnis: - 3x8 Matrix - Aenderungen/Korrekturen: - 11/1995 hh Erste Version - - **************************************************************************/ - -int MGradPhi3D(double* vf, double r, double s, double t) -{ - return MGradOmega3D(vf, r, s, t); -} - -/*************************************************************************** - ROCKFLOW - Funktion: MGetCoor - Aufgabe: - lokale Koordinaten der Eckpunkte - Formalparameter: - E: typ Element-Dimension - 1 - E: nr lokale Knotennummer - X: r,s,t lokale Koordinaten - Ergebnis: - -void- - Aenderungen/Korrekturen: - 04/1996 cb Erste Version - - **************************************************************************/ -void MGetCoor(int typ, long j, double* r, double* s, double* t) -{ - switch (typ) - { - case 0: - *r = (double)(1l - j - j); - *s = *t = 0.0; - break; - case 1: - *r = 1.0 - 2.0 * (double)(((j + 1l) % 4) / 2l); - *s = 1.0 - 2.0 * (double)(j / 2l); - *t = 0.0; - break; - case 2: - *r = 1.0 - 2.0 * (double)(((j + 1l) % 4) / 2l); - *s = 1.0 - 2.0 * (double)((j % 4) / 2l); - *t = 1.0 - 2.0 * (double)(j / 4l); - break; - } /* switch typ */ -} - -/************************************************************************** - ROCKFLOW - Function: MAngleVectors - - Task: Calculate angle between 2 vectors - - Parameter: (I: Input; R: Return; X: Both) - I: *v1, *v2 - - Return: - Angle - - Programming: - 09/2002 MB First Version -**************************************************************************/ -double MAngleVectors(double* v1, double* v2) -{ - double Pproduct; - double angle; - - Pproduct = v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]; - angle = (acos(Pproduct / (MBtrgVec(v1, 3) * MBtrgVec(v2, 3)))) * 180 / PI; - - return angle; -} - -/*************************************************************************** - ROCKFLOW - Funktion: MNormiere - Aufgabe: - Berechnet Norm von Vektor - Formalparameter: - E: *vec - E: n - Ergebnis: - normierter Vektor - Aenderungen/Korrekturen: - 07/1995 hh Erste Version - 11/1999 C.Thorenz Register-Variablen - **************************************************************************/ - -void MNormiere(double* vec, long n) -{ - register long i; - register double vl; - vl = MBtrgVec(vec, n); - for (i = 0; i < n; i++) - vec[i] = vec[i] / vl; -} - -/*************************************************************************** - ROCKFLOW - Funktion: M2Determinante - Aufgabe: - Berechnung der Determinante einer 2x2-Matrix - Formalparameter: - E: *matrix - Ergebnis: - Determinante - Aenderungen/Korrekturen: - 07/1994 hh Erste Version - - **************************************************************************/ - -double M2Determinante(double* matrix) -{ - return matrix[0] * matrix[3] - matrix[1] * matrix[2]; -} /* extern double M2Determinante */ - -/*************************************************************************** - ROCKFLOW - Funktion: Mxg2Determinante - Aufgabe: - Berechnung der Determinante einer mxn-Matrix (n=m und n>2) - --- n --- ---- n --------- - |********* | 0 1 2 3 4 ... n - |********* oder | n+1 n+2 ... - m********* m - |********* | ... - |********* | ... n*m-1 - Formalparameter: - E: *matrix - E: m,n - Ausdehnung - Ergebnis: - Determinante oder 0.0 bei Fehler - Aenderungen/Korrekturen: - 07/1994 hh Erste Version - 11/1999 C.Thorenz Register-Variablen - **************************************************************************/ - -double Mxg2Determinante(double* matrix, long m, long n) -{ - register long i, k, sprung; - register double depp = 0.0, dussel; -#ifdef ERROR_CONTROL - if (m != n) - { - DisplayErrorMsg("Determinate einer nicht-quadratischen Matrix ?"); - return 0.0; - } /* if */ - if (m < 3) - { - DisplayErrorMsg("Determinate zu klein m>2 !!"); - return 0.0; - } /* if */ -#endif - - dussel = 1.0; - for (i = 0; i < m * n; i += m + 1) - dussel *= matrix[i]; - depp += dussel; - for (k = 1; k < m; k++) - { - dussel = 1.0; - sprung = (m - k) * n - 1; - for (i = k; i < n * m; i += m + 1) - { - dussel *= matrix[i]; - if (sprung == i) - i -= m; - } /* for */ - depp += dussel; - } /* for */ - dussel = 1.0; - for (i = n - 1; i < m * n - m + 1; i += m - 1) - dussel *= matrix[i]; - depp -= dussel; - for (k = n - 2; k > -1; k--) - { - dussel = 1.0; - sprung = k * n; - for (i = k; i < n * m; i += m - 1) - { - dussel *= matrix[i]; - if (sprung == i) - i += m; - } /* for */ - depp -= dussel; - } /* for */ - return depp; -} /* extern double Mxg2Determinante */ - -/************************************************************************** - ROCKFLOW - Funktion: MTranspoVec - Aufgabe: - Vektor transponieren - Formalparameter: - X: *vec - E: g - Ergebnis: - Transponierte - Aenderungen/Korrekturen: - 09/1994 hh Erste Version - 11/1999 C.Thorenz Register-Variablen -**************************************************************************/ - -void MTranspoVec(double* vec, long g) -{ - register long i; - register double zwiebel; - for (i = 0; i < g / 2; i++) - { - zwiebel = vec[i]; - vec[i] = vec[g - 1 - i]; - vec[g - 1 - i] = zwiebel; - } -} /* MTranspoVec */ - -/************************************************************************** - ROCKFLOW - Funktion: MTranspoMat - Aufgabe: - Matrizen transponieren - Formalparameter: - E: *mat1 - E: m=Zeilenzahl - E: n=Spaltenzahl - X: *mat2 - Ergebnis: - Transponierte - Aenderungen/Korrekturen: - 08/1994 hh Erste Version - 08/1995 cb (so geht das irgendwie nicht) - 11/1999 C.Thorenz Register-Variablen -**************************************************************************/ - -void MTranspoMat(double* mat1, long m, long n, double* mat2) -{ - register long i, k; - for (i = 0; i < n; i++) - for (k = 0; k < m; k++) - mat2[i * m + k] = mat1[k * n + i]; -} /* MTranspoMat */ - -/************************************************************************** - ROCKFLOW - Funktion: M2InvertiereUndTransponiere - Aufgabe: - 2x2 Matrizen invertieren |a b| - |c d| - 1 | d -b | - inv m = ---------- | | - det(m) | -c a | - Formalparameter: - X: *m - Ergebnis: - Invertierte und Transponierte - Aenderungen/Korrekturen: - 08/1994 hh Erste Version - 11/1999 C.Thorenz Register-Variablen - 08/2001 MK M2Invertiere M2InvertiereUndTransponiere umbenannt -**************************************************************************/ - -void M2InvertiereUndTransponiere(double* m) -{ - register double eddet, zecke; - register int i; - eddet = m[0] * m[3] - m[1] * m[2]; - if (fabs(eddet) > MKleinsteZahl) - eddet = 1.0 / eddet; - zecke = m[0]; - m[0] = m[3]; - m[3] = zecke; - zecke = -m[1]; - m[1] = -m[2]; /* hier wird auch transponiert */ - m[2] = zecke; - for (i = 0; i < 4; i++) - m[i] *= eddet; -} /* M2InvertiereUndTransponiere */ - -/************************************************************************** - ROCKFLOW - Funktion: M2Invertiere - Aufgabe: - 2x2 Matrizen invertieren |a b| - |c d| - 1 | d -b | - inv m = ---------- | | - det(m) | -c a | - Formalparameter: - X: *m - Ergebnis: - Invertierte - Aenderungen/Korrekturen: - 08/1994 hh Erste Version - 11/1999 C.Thorenz Register-Variablen - 08/2001 M. Kohlmeier Funktion transponierte und invertierte jetzt - invertiert sie nur. Alle Aufrufe ersetzt durch - M2InvertiereUndTransponiere - -**************************************************************************/ - -void M2Invertiere(double* m) -{ - register double eddet, zecke; - register int i; - eddet = m[0] * m[3] - m[1] * m[2]; - if (fabs(eddet) > MKleinsteZahl) - eddet = 1.0 / eddet; - zecke = m[0]; - m[0] = m[3]; - m[3] = zecke; - m[1] = -m[1]; - m[2] = -m[2]; - for (i = 0; i < 4; i++) - m[i] *= eddet; -} /* M2Invertiere */ - -/************************************************************************** - ROCKFLOW - Funktion: M3Invertiere - Aufgabe: - 3x3 Matrizen invertieren - - Formalparameter: - X: *m - Ergebnis: - Invertierte - Aenderungen/Korrekturen: - 06/1996 cb Erste Version nach RaRa - 11/1999 C.Thorenz Register-Variablen -**************************************************************************/ - -void M3Invertiere(double* m) -{ - double z[9]; - register double d; - register int i; - d = m[0] * (m[4] * m[8] - m[7] * m[5]) + m[3] * (m[7] * m[2] - m[1] * m[8]) + m[6] * (m[1] * m[5] - m[4] * m[2]); - if (fabs(d) > MKleinsteZahl) - d = 1.0 / d; - z[0] = m[4] * m[8] - m[7] * m[5]; - z[3] = m[6] * m[5] - m[3] * m[8]; - z[6] = m[3] * m[7] - m[6] * m[4]; - z[1] = m[7] * m[2] - m[1] * m[8]; - z[4] = m[0] * m[8] - m[6] * m[2]; - z[7] = m[1] * m[6] - m[0] * m[7]; - z[2] = m[1] * m[5] - m[4] * m[2]; - z[5] = m[3] * m[2] - m[0] * m[5]; - z[8] = m[0] * m[4] - m[3] * m[1]; - for (i = 0; i < 9; i++) - m[i] = d * z[i]; -} /* M3Invertiere */ - -/************************************************************************** - ROCKFLOW - Funktion: MInvertiere - Aufgabe: - Matrizen invertieren (Hauptdiagonalelemente != Null) - (Verfahren aus Schneider) - Formalparameter: - X: *matrix - E: m,n - Ergebnis: - Invertierte - Aenderungen/Korrekturen: - 04/1996 cb Erste Version - 24.06.1999 OK Warnung bei nichtquadratischen Matrizen - 11/1999 C.Thorenz Register-Variablen, Warnung in #ifdef - -**************************************************************************/ -void MInvertiere(double* mat, long m, long n) -{ - register long i, j, k; - -#ifdef ERROR_CONTROL - if (m != n) - { - DisplayErrorMsg("Abbruch !!! Nicht-quadratische Matrix *** MInvertiere"); - abort(); - } - for (k = 0; k < n; k++) - if (fabs(mat[n * k + k]) < MKleinsteZahl) - { - DisplayErrorMsg("Abbruch !!! Diagonalelement Null in MInvertiere"); - abort(); - } -#endif - m = n; /* ah */ - for (k = 0; k < n; k++) - { - for (j = 0; j < n; j++) - if (j != k) - { - mat[n * k + j] /= (-mat[n * k + k]); - for (i = 0; i < n; i++) - if (i != k) - mat[n * i + j] += (mat[n * k + j] * mat[n * i + k]); - } - mat[n * k + k] = 1.0 / mat[n * k + k]; - for (i = 0; i < n; i++) - if (i != k) - mat[n * i + k] *= mat[n * k + k]; - } -} /* MInvertiere */ - -/************************************************************************** - ROCKFLOW - Funktion: MAddVektoren - Aufgabe: - Vektorenen addieren - Formalparameter: - E: *v1, *v2 - A: *vout - E: g - Ergebnis: - Vektorsumme - Aenderungen/Korrekturen: - 08/1994 hh Erste Version - 11/1999 C.Thorenz Register-Variablen - 10/2005 PA Vectorization -**************************************************************************/ -#ifdef SX -int MAddVektoren(double* restrict v1, double* restrict v2, double* restrict vout, long g) -#else -int MAddVektoren(double* v1, double* v2, double* vout, long g) -#endif -{ - register long i; -#ifdef SX -#pragma cdir nodep -#endif - for (i = 0; i < g; i++) - vout[i] = v1[i] + v2[i]; - return 1; -} /* MAddVektoren */ - -/************************************************************************** - ROCKFLOW - Funktion: MAddMatrizen - Aufgabe: - Matrizen addieren - Formalparameter: - E: *m1, *m2 - A: *mout - E: m,n - Ergebnis: - Matrixsumme - Aenderungen/Korrekturen: - 08/1994 hh Erste Version - 11/1999 C.Thorenz Register-Variablen -**************************************************************************/ -#ifdef SX -int MAddMatrizen(double* restrict m1, double* restrict m2, double* restrict mout, long m, long n) -#else -int MAddMatrizen(double* m1, double* m2, double* mout, long m, long n) -#endif -{ - register long i; -#ifdef SX -#pragma cdir nodep -#endif - for (i = 0; i < m * n; i++) - mout[i] = m1[i] + m2[i]; - return 1; -} - -/************************************************************************** - ROCKFLOW - Funktion: MMultVecSkalar - Aufgabe: - Multiplikation eines Vektors mit einem Skalarwert, - wobei der Vektor (!) veraendert wird - Formalparameter: - X: *vec - Zeiger auf Vektor - E: skal - Skalarwert - E: g - Ausdehnung - Aenderungen/Korrekturen: - 07/1994 hh Erste Version - 11/1999 C.Thorenz Register-Variablen -**************************************************************************/ -#ifdef SX -int MMultVecSkalar(double* restrict vec, double skal, long g) -#else -int MMultVecSkalar(double* vec, double skal, long g) -#endif -{ - register long i; -#ifdef SX -#pragma cdir nodep -#endif - for (i = 0; i < g; i++) - vec[i] *= skal; - return 1; -} /* extern int MMultVecSkalar */ - -/************************************************************************** - ROCKFLOW - Funktion: MMultMatSkalar - Aufgabe: - Multiplikation einer Matrix mit einem Skalarwert, - wobei die Matrix (!) veraendert wird - Formalparameter: - X: *matrix - Zeiger auf Matrix - E: skal - Skalarwert - E: m,n - Ausdehnung - Aenderungen/Korrekturen: - 07/1994 hh Erste Version - 11/1999 C.Thorenz Register-Variablen -**************************************************************************/ -#ifdef SX -int MMultMatSkalar(double* restrict matrix, double skal, long m, long n) -#else -int MMultMatSkalar(double* matrix, double skal, long m, long n) -#endif -{ - register long i; -#ifdef SX -#pragma cdir nodep -#endif - for (i = 0; i < m * n; i++) - matrix[i] *= skal; - return 1; -} - -/*########################################################################## - Auf Matrizen und Vektoren aufbauende Funktionen - ######################################################################## */ - -/************************************************************************** - ROCKFLOW - Funktion: MKTF2Dr2D - Aufgabe: - Koordinaten-Transformation 2D -> 2D (r=rotiert) - Formalparameter: - E: *vec1 - E: winkel - A: *vec2 - Ergebnis: - gedrehter Vektor vec2 - Transformation von System 1 zu System 2 - x1 - /\ / x2 -| / um winkel verdrehtes Koordinatensystem -| / es wird immer nach rechts gedreht -| / -|||/ - -+-----------> -|||.. y1 - ... y2 - - Aenderungen/Korrekturen: - 07/1995 hh Erste Version -**************************************************************************/ - -int MKTF2Dr2D(double* vec1, double winkel, double* vec2) -{ -#ifdef ERROR_CONTROL - if (vec2 == NULL) - { - DisplayErrorMsg("MKTF2Dr2D:vec2 muss mit malloc/MMachMat def. sein"); - return 0; - } - if (vec1 == NULL) - { - DisplayErrorMsg("MKTF2Dr2D:vec1 muss mit malloc/MMachMat def. sein"); - return 0; - } -#endif - vec2[0] = vec1[0] * cos(winkel) + vec1[1] * sin(winkel); - vec2[1] = -vec1[0] * sin(winkel) + vec1[1] * cos(winkel); - return 1; -} - -/************************************************************************** - ROCKFLOW - Funktion: MKTFMat2Dr2D - Aufgabe: - Berchnung einer Transformationsmatrix fuer die - Koordinaten-Transformation 2D -> 2D (r=rotiert) - Formalparameter: - E: winkel - A: *tmat - Ergebnis: - Transformationmatrix von System 1 zu System 2 - x1 - /\ / x2 -| / um winkel verdrehtes Koordinatensystem -| / es wird immer nach rechts gedreht -| / -|||/ - -+-----------> -|||.. y1 - ... y2 - - Aenderungen/Korrekturen: - 07/1995 hh Erste Version -**************************************************************************/ - -int MKTFMat2Dr2D(double winkel, double* tmat) -{ -#ifdef ERROR_CONTROL - if (tmat == NULL) - { - DisplayErrorMsg("MKTFMat2Dr2D:tmat muss mit malloc/MMachMat def. sein"); - return 0; - } -#endif - tmat[0] = cos(winkel); - /* tmat[1]= sin(winkel); - tmat[2]=-sin(winkel); */ - tmat[1] = -sin(winkel); - tmat[2] = sin(winkel); - tmat[3] = cos(winkel); - return 1; -} - -/************************************************************************** - ROCKFLOW - Funktion: MKTF2Dt2D - Aufgabe: - Koordinaten-Transformation 2D -> 2D (t=transversal) - Formalparameter: - E: *vec1 - E: dx,dy - A: *vec2 - Ergebnis: - verschobener Vektor vec2 - Transformation von System 1 zu System 2 - x1 - /\ /\ x2 -| | -|||<--->| -| dx | -| +-------------> - -+-----|-----> y2 -| | y1 - - Aenderungen/Korrekturen: - 07/1995 hh Erste Version -**************************************************************************/ - -int MKTF2Dt2D(double* vec1, double dx, double dy, double* vec2) -{ - vec2[1] = vec1[1] + dx; - vec2[2] = vec1[2] + dy; - return 1; -} - -/************************************************************************** - ROCKFLOW - Funktion: MKTF3Dr2D - Aufgabe: - Koordinaten-Transformation 3D-Ebene -> 2D (r=rotiert) - Formalparameter: - E: *vec1 - E: *vec2 - E: winkel - A: *vec - Ergebnis: - gedrehter Vektor vec - Transformation von System 1 (in 3D-Ebene) - zu System 2 (nur noch 2D) - x1 - /\ / x2 -| / um winkel verdrehtes Koordinatensystem -| / es wird immer nach rechts gedreht -| / -|||/ - -+-----------> -|||.. y1 - ... y2 - - Augegeben wird der gedrehte Vektor vec1 - - Aenderungen/Korrekturen: - 07/1995 hh Erste Version - 11/1995 msr Malloc raus / static rein -**************************************************************************/ - -int MKTF3Dr2D(double* vec1, double* vec2, double winkel, double* vec) -{ - static double f; - static double zwick[6]; - /* static double zwack[6]; */ - static double zwock[4]; - static long i; - MNormiere(vec1, 3); - f = MSkalarprodukt(vec1, vec2, 3); - for (i = 0; i < 3; i++) - vec2[i] -= f * vec1[i]; - MNormiere(vec2, 3); - zwick[0] = vec1[0]; - zwick[1] = vec2[0]; - zwick[2] = vec1[1]; - zwick[3] = vec2[1]; - zwick[4] = vec1[2]; - zwick[5] = vec2[2]; - MKTFMat2Dr2D(winkel, zwock); - MMultMatMat(zwick, 3, 2, zwock, 2, 2, vec, 3, 2); - /* MMultMatMat(zwick,3,2,zwock,2,3,zwack,3,2); - zwock[0]=vec[0];zwock[1]=vec[1]; - MMultMatVec(zwack,3,2,zwock,2,vec,2); */ - return 1; -} - -/************************************************************************** - ROCKFLOW - Funktion: MKTFMat3Dr2D - Aufgabe: - Koordinaten-Transformation 3D-Ebene -> 2D (r=rotiert) - Formalparameter: - E: *vec1 - E: *vec2 - E: winkel - A: *vec - Ergebnis: - gedrehter Vektor vec - Transformation von System 1 (in 3D-Ebene) - zu System 2 (nur noch 2D) - x1 - /\ / x2 -| / um winkel verdrehtes Koordinatensystem -| / es wird immer nach rechts gedreht -| / -|||/ - -+-----------> -|||.. y1 - ... y2 - - Augegeben wird der gedrehte Vektor vec1 - - Aenderungen/Korrekturen: - 07/1995 hh Erste Version - 11/1995 msr Malloc raus / static rein -**************************************************************************/ - -void MKTFMat3Dr2D(double* vec1, double* vec2, double winkel, double* mat) -{ /* if(winkel!=0.0) */ /* diese Verzweigung spart viel Zeit */ - { - static double f; - static double zwick[6]; - static double zwack[6]; - static double zwock[4]; - static long i; - MNormiere(vec1, 3); - f = MSkalarprodukt(vec1, vec2, 3); - for (i = 0; i < 3; i++) - vec2[i] -= f * vec1[i]; - MNormiere(vec2, 3); - zwick[0] = vec1[0]; - zwick[1] = vec2[0]; - zwick[2] = vec1[1]; - zwick[3] = vec2[1]; - zwick[4] = vec1[2]; - zwick[5] = vec2[2]; - MKTFMat2Dr2D(winkel, zwock); - MMultMatMat(zwick, 3, 2, zwock, 2, 2, zwack, 3, 2); - for (i = 0; i < 6; i++) - mat[i] = zwack[i]; - } /* if */ - - /* else - {MNulleMat(mat,3,2); - mat[0]=1.0;mat[3]=1.0; - } */ - /* else */ -} - -/*########################################################################## - Pruefunktion fuer Matrix - ######################################################################## */ - -/************************************************************************** - ROCKFLOW - Funktion: MBistDuDiagMat - Aufgabe: - Prueft Matrix auf diagonalitaet - Formalparameter: - E: *matrix - E: m,n - Ergebnis: - -1 - keine quadartische Matrix - 0 - keine Diagonalmatrix - 1 - Diagonalmatrix - 3 - Tridiagonalmatrix - Aenderungen/Korrekturen: - 09/1994 hh Erste Version -**************************************************************************/ - -int MBistDuDiagMat(double* matrix, long m, long n) -{ - long i, j, knall = 0, diag = 0; - if (m != n) - return -1; - for (i = 0; i < m; i++) - for (j = 0; j < n; j++) - /*printf(" %ld %ld \n",i,j); */ - if ((fabs(matrix[i + j * n]) > MKleinsteZahl) && ((i - j < -1) || (i - j > 1))) - return 0; - - diag = 3; - knall = 0; - for (i = 0; i < m * n; i += m + 1) - for (j = i + 1; ((j < i + m) && (j < m * n)); j++) - if (fabs(matrix[j]) > MKleinsteZahl) - knall = 1; - if (!knall) - diag = 1; - return diag; -} /* of MBistDuDiagMat */ - -/*########################################################################## - Bearbeitungfunktionen fuer Vektoren - ######################################################################## */ - -/************************************************************************** - ROCKFLOW - Funktion: MLoeschVec - Aufgabe: - Zerstoeren eines bestehenden Vektors - Formalparameter: - E: *vec - Ergebnis: - Rueckgabewert 1 - Aenderungen/Korrekturen: - 08/1994 hh Erste Version -**************************************************************************/ - -void MLoeschVec(double* vec) -{ -#ifdef ERROR_CONTROL - if (vec == NULL) - DisplayErrorMsg("Fehler in MLoeschVec"); -#endif - vec = (double*)Free(vec); -} /* MLoeschVec */ - -/*########################################################################## - Bearbeitungfunktionen fuer Matrizen - ######################################################################## */ - -/************************************************************************** - ROCKFLOW - Funktion: MMachMat - Aufgabe: - Erzeugen einer neuen Matrix - Formalparameter: - E: m,n - Ergebnis: - Zeiger auf die erzeugte Matrix. - oder Rueckgabewert NULL wenn nicht genuegend Speicher da ist. - Aenderungen/Korrekturen: - 08/1994 hh Erste Version - 11/1999 C.Thorenz Nullung herausgenommen -**************************************************************************/ - -double* MMachMat(long m, long n) -{ - double* zwerg; - zwerg = (double*)Malloc(sizeof(double) * m * n); -#ifdef ERROR_CONTROL - if (zwerg == NULL) - DisplayErrorMsg("zuwenig Speicher in MMachMat"); -#endif - return zwerg; -} /* MMachMat */ - -/************************************************************************** - ROCKFLOW - Funktion: MNullMat - Aufgabe: - Fuellt Matrix mit Nullen - Formalparameter: - E: *zwerg,m,n - Ergebnis: - Aenderungen/Korrekturen: - 11/1999 C.Thorenz Erste Version -**************************************************************************/ - -void MNullMat(double* zwerg, long m, long n) -{ - register long i; -#ifdef ERROR_CONTROL - if (zwerg == NULL) - DisplayErrorMsg("Fehler in MNullMat"); -#endif - for (i = 0; i < m * n; i++) - zwerg[i] = 0.0; -} /* MNullMat */ - -/************************************************************************** - ROCKFLOW - Funktion: MLoeschMat - Aufgabe: - Zerstoeren einer bestehenden Matrix - Formalparameter: - E: *vec - Ergebnis: - Rueckgabewert 1 - Aenderungen/Korrekturen: - 08/1994 hh Erste Version -**************************************************************************/ - -void MLoeschMat(double* mat) -{ -#ifdef ERROR_CONTROL - if (mat == NULL) - DisplayErrorMsg("Fehler in MLoeschMat"); -#endif - mat = (double*)Free(mat); -} /* MLoeschMat */ - -/************************************************************************** - ROCKFLOW - Funktion: MKopierMat - Aufgabe: - Kopieren einer Matrix - Formalparameter: - E: *matq - A: *mat - E: m - E: n - Ergebnis: - Rueckgabewert 1 - Aenderungen/Korrekturen: - 08/1994 hh Erste Version -**************************************************************************/ - -void MKopierMat(double* matq, double* matz, long m, long n) -{ - register long i; - for (i = 0; i < m * n; i++) - matz[i] = matq[i]; -} /* MKopierMat */ - -/*########################################################################## - Ausgabefunktionen - ########################################################################### */ -/************************************************************************** - ROCKFLOW - Funktion: MZeigVec - Aufgabe: - Anzeigen eines Vektors - Formalparameter: - E: *vec - E: grad - E: *text - Programmaenderungen: - 08/1994 hh Erste Version -**************************************************************************/ -void MZeigVec(double* vec, long grad, char* text) -{ - long i; - printf("\n%s\n", text); - for (i = 0; i < grad; i++) - printf(" | %e | \n", vec[i]); -} - -/************************************************************************** - ROCKFLOW - Funktion: M2FileVec - Aufgabe: - Anzeigen eines Vektors - Formalparameter: - E: *vec - E: grad - E: *text - Programmaenderungen: - 08/1994 hh Erste Version -**************************************************************************/ -void M2FileVec(double* vec, long grad, char* text) -{ - long i; - FILE* fp; - fp = fopen("fgnuplot.asc", "a"); - fprintf(fp, "\n%s\n", text); - for (i = 0; i < grad; i++) - fprintf(fp, "%e\n", vec[i]); - fclose(fp); -} - -/************************************************************************** - ROCKFLOW - Funktion: MZeigMat - Aufgabe: - Anzeigen einer Matrix - Formalparameter: - E: *mat - E: m=Zeilenzahl - E: n=Spaltenzahl - E: *text - Programmaenderungen: - 08/1994 hh Erste Version -**************************************************************************/ -void MZeigMat(double* mat, long m, long n, char* text) -{ - long i, k; - printf("\n%s\n", text); - for (k = 0; k < m; k++) - { - printf("| "); - for (i = 0; i < n; i++) - printf(" %e ", mat[i + k * n]); - puts(" |"); - } -} /* MZeigMat */ - -/************************************************************************** - ROCKFLOW - Funktion: MVekNormMax - - Aufgabe: - Berechnet die Maximumnorm von x - - Formalparameter: (E: Eingabe; R: Rueckgabe; X: Beides) - E double *x : Zeiger auf Vektor - E long n : Dimension von x - - Ergebnis: - Norm - - Programmaenderungen: - 11/1995 MSR Erste Version - -**************************************************************************/ -double MVekNormMax(double* x, long n) -{ - register long i; - register double erg = fabs(x[0]); - for (i = 1l; i < n; i++) - if (fabs(x[i]) > erg) - erg = fabs(x[i]); - return erg; -} - -/*########################################################################## - Geometrische Funktionen im R3 - ######################################################################## */ - -/*************************************************************************** - ROCKFLOW - Funktion: MCalcDistancePointToPoint - - Aufgabe: - Abstand zweier Punkte im R3 - - Formalparameter: - E: *pt1 - Punktkoordinaten - E: *pt2 - Punktkoordinaten - - Ergebnis: - Abstand - - Aenderungen/Korrekturen: - 02/2000 C. Thorenz Erste Version - - **************************************************************************/ - -/*************************************************************************** - ROCKFLOW - Funktion: MCalcProjectionOfPointOnLine - - Aufgabe: - Ermittelt den Projektionspunkt eines Punktes auf eine Linie (Lotpunkt) - - Formalparameter: - E: *pt - Punktkoordinaten - E: *l1 - Punktkoordinaten fuer ersten Punkt auf Linie - E: *l2 - Punktkoordinaten fuer zweiten Punkt auf Linie - R: *proj - Punktkoordinaten fuer Projektionspunkt auf Linie - Ergebnis: - Abstand Punkt-Linie - - Aenderungen/Korrekturen: - 02/2000 C. Thorenz Erste Version - - **************************************************************************/ -double MCalcProjectionOfPointOnLine(double* pt, double* l1, double* l2, double* proj) -{ - int i; - double vec1[3], vec2[3], l, projektionslaenge, abstand; - - /* Vektor Linienpunkt zu Punkt */ - for (i = 0; i < 3; i++) - vec1[i] = pt[i] - l1[i]; - - /* Vektor Linienpunkt zu Linienpunkt */ - for (i = 0; i < 3; i++) - vec2[i] = l2[i] - l1[i]; - - /* Normieren */ - l = MBtrgVec(vec2, 3); -#ifdef ERROR_CONTROL - if (l < MKleinsteZahl) - { - printf("\n Fehler in MCalcProjectionOfPointOnLine: Laenge ist Null !!!"); - exit(1); - } -#endif - for (i = 0; i < 3; i++) - vec2[i] /= (l + MKleinsteZahl); - - projektionslaenge = MSkalarprodukt(vec1, vec2, 3); - - /* Punkt bestimmen */ - for (i = 0; i < 3; i++) - proj[i] = l1[i] + projektionslaenge * vec2[i]; - - abstand = MCalcDistancePointToPoint(proj, pt); - - return abstand; -} - -/*************************************************************************** - ROCKFLOW - Funktion: MCalcProjectionOfPointOnPlane - - Aufgabe: - Ermittelt den Projektionspunkt eines Punktes auf einer Flaeche (Lotpunkt) - - Formalparameter: - E: *pt - Punktkoordinaten - E: *e1 - Punktkoordinaten fuer ersten Punkt auf Ebene - E: *e2 - Punktkoordinaten fuer zweiten Punkt auf Ebene - E: *e3 - Punktkoordinaten fuer dritten Punkt auf Ebene - R: *proj - Punktkoordinaten fuer Projektionspunkt auf Linie - Ergebnis: - Abstand-Flaeche (Kann negativ sein, je nach Lage des Punkts zur Ebene) - - Aenderungen/Korrekturen: - 02/2000 C. Thorenz Erste Version - - **************************************************************************/ - -/*************************************************************************** - ROCKFLOW - Funktion: IsNodeInsideTriangle - - Aufgabe: Ueberprueft, ob sich ein vorgegebener Knoten innerhalb - eines vorgegebenen Dreiecks befindet - - Formalparameter: - E long n1, n2, n3: Knoten des Dreiecks - E long node: Knoten - - Ergebnis: - 0 bei Fehler, sonst 1 - - Aenderungen/Korrekturen: - 05/2000 RK Erste Version - - **************************************************************************/ - -/*************************************************************************** - ROCKFLOW - Funktion: CalcTriangleArea - - Aufgabe: Berechnet die Flaeche eines vorgegebenen Dreiecks befindet - - Formalparameter: - E long n1, n2, n3: Knoten des Dreiecks - - Ergebnis: Dreiecksflaeche - - Aenderungen/Korrekturen: - 05/2000 RK Erste Version - - **************************************************************************/ -double CalcTriangleArea(long n1, long n2, long n3) -{ - double area = 0.0; - n1 = n1; - n2 = n2; - n3 = n3; // OK411 - /*OK411 - double vec1[3], vec2[3], vec3[3]; - - vec1[0] = GetNodeX(n2)-GetNodeX(n1); - vec1[1] = GetNodeY(n2)-GetNodeY(n1); - vec1[2] = GetNodeZ(n2)-GetNodeZ(n1); - - vec2[0] = GetNodeX(n1)-GetNodeX(n3); - vec2[1] = GetNodeY(n1)-GetNodeY(n3); - vec2[2] = GetNodeZ(n1)-GetNodeZ(n3); - - M3KreuzProdukt(vec1,vec2,vec3); - - area = 0.5 * MBtrgVec(vec3,3); - */ - return area; -} - -/*************************************************************************** - ROCKFLOW - Funktion: IsNodeInsidePlain - - Aufgabe: Ueberprueft, ob sich ein vorgegebener Knoten innerhalb - einer vorgegebenen Flaeche befindet - - Formalparameter: - E long n1, n2, n3: Flaeche aufspannenden Knoten - E long node: Knoten - - Ergebnis: - 0 bei Fehler, sonst 1 - - Aenderungen/Korrekturen: - 05/2000 RK Erste Version - - **************************************************************************/ - -/*########################################################################## - Funktionen fuer zweidimensionale 9-Knoten-Elemente - ######################################################################## */ - -/*************************************************************************** - ROCKFLOW - Funktion: MOmega2D_9N - Aufgabe: - Berechnet Ansatzfunktion (r,s). - - Formalparameter: - Z: *vf - 1x9 Feld - E: r,s (s.o.) - Ergebnis: - Vektor - Aenderungen/Korrekturen: - 12/1999 RK Erste Version - - **************************************************************************/ - -int MOmega2D_9N(double* vf, double r, double s) -{ - int ok = 0; -#ifdef ERROR_CONTROL - if (vf == NULL) - { - ok = 0; - return ok; - } -#endif - vf[8] = (1.0 - r * r) * (1.0 - s * s); - vf[7] = 0.5 * (1.0 - s * s) * (1.0 + r) - 0.5 * vf[8]; - vf[6] = 0.5 * (1.0 - r * r) * (1.0 - s) - 0.5 * vf[8]; - vf[5] = 0.5 * (1.0 - s * s) * (1.0 - r) - 0.5 * vf[8]; - vf[4] = 0.5 * (1.0 - r * r) * (1.0 + s) - 0.5 * vf[8]; - vf[3] = 0.25 * (1.0 + r) * (1.0 - s) - 0.5 * vf[6] - 0.5 * vf[7] - 0.25 * vf[8]; - vf[2] = 0.25 * (1.0 - r) * (1.0 - s) - 0.5 * vf[5] - 0.5 * vf[6] - 0.25 * vf[8]; - vf[1] = 0.25 * (1.0 - r) * (1.0 + s) - 0.5 * vf[4] - 0.5 * vf[5] - 0.25 * vf[8]; - vf[0] = 0.25 * (1.0 + r) * (1.0 + s) - 0.5 * vf[4] - 0.5 * vf[7] - 0.25 * vf[8]; - return ok = 1; -} - -/*************************************************************************** - ROCKFLOW - Funktion: MPhi2D_9N - - Aufgabe: Berechnet Testfunktion (r,s) - - Formalparameter: Z: *vf - 1x9 Feld - E: r,s (s.o.) - - Ergebnis: Vektor - - Aenderungen/Korrekturen: - 12/1999 R.Kaiser Erste Version - - **************************************************************************/ - -int MPhi2D_9N(double* vf, double r, double s) -{ - int ok = 0; - - vf[8] = (1.0 - r * r) * (1.0 - s * s); - vf[7] = 0.5 * (1.0 - s * s) * (1.0 + r) - 0.5 * vf[8]; - vf[6] = 0.5 * (1.0 - r * r) * (1.0 - s) - 0.5 * vf[8]; - vf[5] = 0.5 * (1.0 - s * s) * (1.0 - r) - 0.5 * vf[8]; - vf[4] = 0.5 * (1.0 - r * r) * (1.0 + s) - 0.5 * vf[8]; - vf[3] = 0.25 * (1.0 + r) * (1.0 - s) - 0.5 * vf[6] - 0.5 * vf[7] - 0.25 * vf[8]; - vf[2] = 0.25 * (1.0 - r) * (1.0 - s) - 0.5 * vf[5] - 0.5 * vf[6] - 0.25 * vf[8]; - vf[1] = 0.25 * (1.0 - r) * (1.0 + s) - 0.5 * vf[4] - 0.5 * vf[5] - 0.25 * vf[8]; - vf[0] = 0.25 * (1.0 + r) * (1.0 + s) - 0.5 * vf[4] - 0.5 * vf[7] - 0.25 * vf[8]; - return ok = 1; -} - -/*************************************************************************** - ROCKFLOW - Funktion: MGradOmega2D_9N - Aufgabe: - Berechnet Gradient eines 2D Vektorfeldes, - dessen Ansatzfunktion (Vektor) bekannt ist (r,s). - - Formalparameter: - Z: *vf - 2x9 Feld - E: r,s (s.o.) - Ergebnis: - 2x9 Matrix - Aenderungen/Korrekturen: - 12/1999 RK Erste Version - - **************************************************************************/ - -int MGradOmega2D_9N(double* vf, double r, double s) -{ - int ok = 0; -#ifdef ERROR_CONTROL - if (vf == NULL) - { - ok = 0; - return ok; - } -#endif - - vf[8] = -2.0 * r * (1.0 - s * s); - vf[7] = +0.5 * (1.0 - s * s) - 0.5 * vf[8]; - vf[6] = -1.0 * r * (1.0 - s) - 0.5 * vf[8]; - vf[5] = -0.5 * (1.0 - s * s) - 0.5 * vf[8]; - vf[4] = -1.0 * r * (1.0 + s) - 0.5 * vf[8]; - vf[3] = +0.25 * (1 - s) - 0.5 * vf[6] - 0.5 * vf[7] - 0.25 * vf[8]; - vf[2] = -0.25 * (1 - s) - 0.5 * vf[5] - 0.5 * vf[6] - 0.25 * vf[8]; - vf[1] = -0.25 * (1 + s) - 0.5 * vf[4] - 0.5 * vf[5] - 0.25 * vf[8]; - vf[0] = +0.25 * (1 + s) - 0.5 * vf[4] - 0.5 * vf[7] - 0.25 * vf[8]; - - vf[17] = -2.0 * s * (1.0 - r * r); - vf[16] = -1.0 * s * (1.0 + r) - 0.5 * vf[17]; - vf[15] = -0.5 * (1.0 - r * r) - 0.5 * vf[17]; - vf[14] = -1.0 * s * (1.0 - r) - 0.5 * vf[17]; - vf[13] = +0.5 * (1 - r * r) - 0.5 * vf[17]; - vf[12] = -0.25 * (1 + r) - 0.5 * vf[15] - 0.5 * vf[16] - 0.25 * vf[17]; - vf[11] = -0.25 * (1 - r) - 0.5 * vf[14] - 0.5 * vf[15] - 0.25 * vf[17]; - vf[10] = +0.25 * (1 - r) - 0.5 * vf[13] - 0.5 * vf[14] - 0.25 * vf[17]; - vf[9] = +0.25 * (1 + r) - 0.5 * vf[13] - 0.5 * vf[16] - 0.25 * vf[17]; - - return ok = 1; -} - -/*************************************************************************** - ROCKFLOW - Funktion: MGradPhi2D_9N - Aufgabe: - Berechnet Gradient der Testfunktionen (r,s) - grad phi = grad omega - Formalparameter: - Z: *vf - 2x9 Feld - E: r,s (s.o.) - Ergebnis: - 2x9 Matrix - Aenderungen/Korrekturen: - 12/1999 RK Erste Version - - **************************************************************************/ - -int MGradPhi2D_9N(double* vf, double r, double s) -{ - return MGradOmega2D_9N(vf, r, s); -} - -/*########################################################################## - Funktionen fuer dreidimensionale 20-Knoten-Elemente - ######################################################################## */ - -/*************************************************************************** - ROCKFLOW - Funktion: MPhi3D_20N - Aufgabe: - Berechnet Testfunktion (r,s,t). - (K.- J. Bathe, Finite-Elemente-Methoden) - - Formalparameter: - Z: *vf - 1x20 Feld - E: r,s,t (s.o.) - Ergebnis: - Vektor - Aenderungen/Korrekturen: - 10/2001 RK Erste Version - - **************************************************************************/ -int MPhi3D_20N(double* vf, double r, double s, double t) -{ - int ok = 0; - double g[20]; - -#ifdef ERROR_CONTROL - if (vf == NULL) - { - ok = 0; - return ok; - } -#endif - - g[0] = 0.125 * (1.0 + r) * (1.0 + s) * (1.0 + t); - g[1] = 0.125 * (1.0 - r) * (1.0 + s) * (1.0 + t); - g[2] = 0.125 * (1.0 - r) * (1.0 - s) * (1.0 + t); - g[3] = 0.125 * (1.0 + r) * (1.0 - s) * (1.0 + t); - g[4] = 0.125 * (1.0 + r) * (1.0 + s) * (1.0 - t); - g[5] = 0.125 * (1.0 - r) * (1.0 + s) * (1.0 - t); - g[6] = 0.125 * (1.0 - r) * (1.0 - s) * (1.0 - t); - g[7] = 0.125 * (1.0 + r) * (1.0 - s) * (1.0 - t); - g[8] = 0.25 * (1.0 - r * r) * (1.0 + s) * (1.0 + t); - g[9] = 0.25 * (1.0 - r) * (1.0 - s * s) * (1.0 + t); - g[10] = 0.25 * (1.0 - r * r) * (1.0 - s) * (1.0 + t); - g[11] = 0.25 * (1.0 + r) * (1.0 - s * s) * (1.0 + t); - g[12] = 0.25 * (1.0 - r * r) * (1.0 + s) * (1.0 - t); - g[13] = 0.25 * (1.0 - r) * (1.0 - s * s) * (1.0 - t); - g[14] = 0.25 * (1.0 - r * r) * (1.0 - s) * (1.0 - t); - g[15] = 0.25 * (1.0 + r) * (1.0 - s * s) * (1.0 - t); - g[16] = 0.25 * (1.0 + r) * (1.0 + s) * (1.0 - t * t); - g[17] = 0.25 * (1.0 - r) * (1.0 + s) * (1.0 - t * t); - g[18] = 0.25 * (1.0 - r) * (1.0 - s) * (1.0 - t * t); - g[19] = 0.25 * (1.0 + r) * (1.0 - s) * (1.0 - t * t); - - vf[0] = g[0] - (g[8] + g[11] + g[16]) * 0.5; - vf[1] = g[1] - (g[8] + g[9] + g[17]) * 0.5; - vf[2] = g[2] - (g[9] + g[10] + g[18]) * 0.5; - vf[3] = g[3] - (g[10] + g[11] + g[19]) * 0.5; - vf[4] = g[4] - (g[12] + g[15] + g[16]) * 0.5; - vf[5] = g[5] - (g[12] + g[13] + g[17]) * 0.5; - vf[6] = g[6] - (g[13] + g[14] + g[18]) * 0.5; - vf[7] = g[7] - (g[14] + g[15] + g[19]) * 0.5; - - vf[8] = g[8]; - vf[9] = g[9]; - vf[10] = g[10]; - vf[11] = g[11]; - vf[12] = g[12]; - vf[13] = g[13]; - vf[14] = g[14]; - vf[15] = g[15]; - vf[16] = g[16]; - vf[17] = g[17]; - vf[18] = g[18]; - vf[19] = g[19]; - - return ok = 1; -} - -/*************************************************************************** - ROCKFLOW - Funktion: MGradOmega3D_20N - Aufgabe: - Berechnet Gradient der Ansatzfunktion (r,s,t) - - Formalparameter: - Z: *vf - 3x20 Feld - E: r,s,t (s.o.) - Ergebnis: - Vektor - Aenderungen/Korrekturen: - 10/2001 RK Erste Version - - **************************************************************************/ -int MGradOmega3D_20N(double* vf, double r, double s, double t) -{ - int ok = 0; - double g_r[20]; /* r-Ableitung */ - double g_s[20]; /* s-Ableitung */ - double g_t[20]; /* t-Ableitung */ - -#ifdef ERROR_CONTROL - if (vf == NULL) - { - ok = 0; - return ok; - } -#endif - - g_r[0] = +0.125 * (1.0 + s) * (1.0 + t); - g_r[1] = -0.125 * (1.0 + s) * (1.0 + t); - g_r[2] = -0.125 * (1.0 - s) * (1.0 + t); - g_r[3] = +0.125 * (1.0 - s) * (1.0 + t); - g_r[4] = +0.125 * (1.0 + s) * (1.0 - t); - g_r[5] = -0.125 * (1.0 + s) * (1.0 - t); - g_r[6] = -0.125 * (1.0 - s) * (1.0 - t); - g_r[7] = +0.125 * (1.0 - s) * (1.0 - t); - g_r[8] = -0.5 * r * (1.0 + s) * (1.0 + t); - g_r[9] = -0.25 * (1.0 - s * s) * (1.0 + t); - g_r[10] = -0.5 * r * (1.0 - s) * (1.0 + t); - g_r[11] = +0.25 * (1.0 - s * s) * (1.0 + t); - g_r[12] = -0.5 * r * (1.0 + s) * (1.0 - t); - g_r[13] = -0.25 * (1.0 - s * s) * (1.0 - t); - g_r[14] = -0.5 * r * (1.0 - s) * (1.0 - t); - g_r[15] = +0.25 * (1.0 - s * s) * (1.0 - t); - g_r[16] = +0.25 * (1.0 + s) * (1.0 - t * t); - g_r[17] = -0.25 * (1.0 + s) * (1.0 - t * t); - g_r[18] = -0.25 * (1.0 - s) * (1.0 - t * t); - g_r[19] = +0.25 * (1.0 - s) * (1.0 - t * t); - - g_s[0] = +0.125 * (1.0 + r) * (1.0 + t); - g_s[1] = +0.125 * (1.0 - r) * (1.0 + t); - g_s[2] = -0.125 * (1.0 - r) * (1.0 + t); - g_s[3] = -0.125 * (1.0 + r) * (1.0 + t); - g_s[4] = +0.125 * (1.0 + r) * (1.0 - t); - g_s[5] = +0.125 * (1.0 - r) * (1.0 - t); - g_s[6] = -0.125 * (1.0 - r) * (1.0 - t); - g_s[7] = -0.125 * (1.0 + r) * (1.0 - t); - g_s[8] = +0.25 * (1.0 - r * r) * (1.0 + t); - g_s[9] = -0.5 * (1.0 - r) * s * (1.0 + t); - g_s[10] = -0.25 * (1.0 - r * r) * (1.0 + t); - g_s[11] = -0.5 * (1.0 + r) * s * (1.0 + t); - g_s[12] = +0.25 * (1.0 - r * r) * (1.0 - t); - g_s[13] = -0.5 * (1.0 - r) * s * (1.0 - t); - g_s[14] = -0.25 * (1.0 - r * r) * (1.0 - t); - g_s[15] = -0.5 * (1.0 + r) * s * (1.0 - t); - g_s[16] = +0.25 * (1.0 + r) * (1.0 - t * t); - g_s[17] = +0.25 * (1.0 - r) * (1.0 - t * t); - g_s[18] = -0.25 * (1.0 - r) * (1.0 - t * t); - g_s[19] = -0.25 * (1.0 + r) * (1.0 - t * t); - - g_t[0] = +0.125 * (1.0 + r) * (1.0 + s); - g_t[1] = +0.125 * (1.0 - r) * (1.0 + s); - g_t[2] = +0.125 * (1.0 - r) * (1.0 - s); - g_t[3] = +0.125 * (1.0 + r) * (1.0 - s); - g_t[4] = -0.125 * (1.0 + r) * (1.0 + s); - g_t[5] = -0.125 * (1.0 - r) * (1.0 + s); - g_t[6] = -0.125 * (1.0 - r) * (1.0 - s); - g_t[7] = -0.125 * (1.0 + r) * (1.0 - s); - g_t[8] = +0.25 * (1.0 - r * r) * (1.0 + s); - g_t[9] = +0.25 * (1.0 - r) * (1.0 - s * s); - g_t[10] = +0.25 * (1.0 - r * r) * (1.0 - s); - g_t[11] = +0.25 * (1.0 + r) * (1.0 - s * s); - g_t[12] = -0.25 * (1.0 - r * r) * (1.0 + s); - g_t[13] = -0.25 * (1.0 - r) * (1.0 - s * s); - g_t[14] = -0.25 * (1.0 - r * r) * (1.0 - s); - g_t[15] = -0.25 * (1.0 + r) * (1.0 - s * s); - g_t[16] = -0.5 * (1.0 + r) * (1.0 + s) * t; - g_t[17] = -0.5 * (1.0 - r) * (1.0 + s) * t; - g_t[18] = -0.5 * (1.0 - r) * (1.0 - s) * t; - g_t[19] = -0.5 * (1.0 + r) * (1.0 - s) * t; - - vf[0] = g_r[0] - (g_r[8] + g_r[11] + g_r[16]) * 0.5; - vf[1] = g_r[1] - (g_r[8] + g_r[9] + g_r[17]) * 0.5; - vf[2] = g_r[2] - (g_r[9] + g_r[10] + g_r[18]) * 0.5; - vf[3] = g_r[3] - (g_r[10] + g_r[11] + g_r[19]) * 0.5; - vf[4] = g_r[4] - (g_r[12] + g_r[15] + g_r[16]) * 0.5; - vf[5] = g_r[5] - (g_r[12] + g_r[13] + g_r[17]) * 0.5; - vf[6] = g_r[6] - (g_r[13] + g_r[14] + g_r[18]) * 0.5; - vf[7] = g_r[7] - (g_r[14] + g_r[15] + g_r[19]) * 0.5; - - vf[8] = g_r[8]; - vf[9] = g_r[9]; - vf[10] = g_r[10]; - vf[11] = g_r[11]; - vf[12] = g_r[12]; - vf[13] = g_r[13]; - vf[14] = g_r[14]; - vf[15] = g_r[15]; - vf[16] = g_r[16]; - vf[17] = g_r[17]; - vf[18] = g_r[18]; - vf[19] = g_r[19]; - - vf[20] = g_s[0] - (g_s[8] + g_s[11] + g_s[16]) * 0.5; - vf[21] = g_s[1] - (g_s[8] + g_s[9] + g_s[17]) * 0.5; - vf[22] = g_s[2] - (g_s[9] + g_s[10] + g_s[18]) * 0.5; - vf[23] = g_s[3] - (g_s[10] + g_s[11] + g_s[19]) * 0.5; - vf[24] = g_s[4] - (g_s[12] + g_s[15] + g_s[16]) * 0.5; - vf[25] = g_s[5] - (g_s[12] + g_s[13] + g_s[17]) * 0.5; - vf[26] = g_s[6] - (g_s[13] + g_s[14] + g_s[18]) * 0.5; - vf[27] = g_s[7] - (g_s[14] + g_s[15] + g_s[19]) * 0.5; - - vf[28] = g_s[8]; - vf[29] = g_s[9]; - vf[30] = g_s[10]; - vf[31] = g_s[11]; - vf[32] = g_s[12]; - vf[33] = g_s[13]; - vf[34] = g_s[14]; - vf[35] = g_s[15]; - vf[36] = g_s[16]; - vf[37] = g_s[17]; - vf[38] = g_s[18]; - vf[39] = g_s[19]; - - vf[40] = g_t[0] - (g_t[8] + g_t[11] + g_t[16]) * 0.5; - vf[41] = g_t[1] - (g_t[8] + g_t[9] + g_t[17]) * 0.5; - vf[42] = g_t[2] - (g_t[9] + g_t[10] + g_t[18]) * 0.5; - vf[43] = g_t[3] - (g_t[10] + g_t[11] + g_t[19]) * 0.5; - vf[44] = g_t[4] - (g_t[12] + g_t[15] + g_t[16]) * 0.5; - vf[45] = g_t[5] - (g_t[12] + g_t[13] + g_t[17]) * 0.5; - vf[46] = g_t[6] - (g_t[13] + g_t[14] + g_t[18]) * 0.5; - vf[47] = g_t[7] - (g_t[14] + g_t[15] + g_t[19]) * 0.5; - - vf[48] = g_t[8]; - vf[49] = g_t[9]; - vf[50] = g_t[10]; - vf[51] = g_t[11]; - vf[52] = g_t[12]; - vf[53] = g_t[13]; - vf[54] = g_t[14]; - vf[55] = g_t[15]; - vf[56] = g_t[16]; - vf[57] = g_t[17]; - vf[58] = g_t[18]; - vf[59] = g_t[19]; - - return ok = 1; -} - -/*************************************************************************** - ROCKFLOW - Funktion: MGradPhi3D_20N - Aufgabe: - Berechnet Gradient der Testfunktionen (r,s,t) - grad phi = grad omega - Formalparameter: - Z: *vf - 3x20 Feld - E: r,s (s.o.) - Ergebnis: - 3x20 Matrix - Aenderungen/Korrekturen: - 10/2001 RK Erste Version - - **************************************************************************/ -int MGradPhi3D_20N(double* vf, double r, double s, double t) -{ - return MGradOmega3D_20N(vf, r, s, t); -} - -/*########################################################################## - Sortierfunktionen - ########################################################################*/ -/*************************************************************************** - ROCKFLOW - Funktion: MCompare_for_MQSort_LongDouble - Aufgabe: - Vergleichfunktion fuer zugehoerige Sortierfunktion - Formalparameter: - *Arg1, *Arg2 - - Ergebnis: - 1,0,-1 (groesser, gleich, kleiner) - Aenderungen/Korrekturen: - 06/2001 MK Erste Version - - **************************************************************************/ -int MCompare_for_MQSort_LongDouble(const void* Arg1, const void* Arg2) -{ - typedef struct - { - long index; - double value; - } MType_LongDouble; - MType_LongDouble *arg1, *arg2; - arg1 = (MType_LongDouble*)Arg1; - arg2 = (MType_LongDouble*)Arg2; - if ((arg1->value) > (arg2->value)) - return 1; - if (fabs(arg1->value - arg2->value) < MKleinsteZahl) - return 0; - else - return -1; -} - -/*************************************************************************** - ROCKFLOW - Funktion: MQSort_LongDouble - Aufgabe: - Sortierfunktion - Formalparameter: - *DataSets - NumberOfDataSets - SizeOfDataSet - - Ergebnis: - Sortierter Datensatz - Aenderungen/Korrekturen: - 06/2001 MK Erste Version - - **************************************************************************/ -void MQSort_LongDouble(void* DataSets, const int NumberOfDataSets, const int SizeOfDataSet) -{ - qsort(DataSets, NumberOfDataSets, SizeOfDataSet, MCompare_for_MQSort_LongDouble); -} - -/*########################################################################## - eventuell noch nuetzliche Funktionen, - ######################################################################## */ - -/*************************************************************************** - ROCKFLOW - Funktion: JWDMMultMatSkalar - - Aufgabe: - Multiplikation einer Matrix mit einem Skalarwert, - wobei die Matrix nicht veraendert wird - Formalparameter: - E: *matrix - Zeiger auf Matrix - E: skal - Skalarwert - E: m,n - Ausdehnung - Ergebnis: - Zeiger auf Ergebnis - Aenderungen/Korrekturen: - 07/1994 hh Erste Version - - **************************************************************************/ - -double* JWDMMultMatSkalar(double* matrix, double skal, long m, long n) -{ /* flip ist Zwischenspeicher fuer die Matrix */ - double *flip, *flippy; - long i; - flip = (double*)Malloc(sizeof(double) * n * m); - flippy = flip; - if (flip == NULL) - { - DisplayErrorMsg(" kein Speicher mehr !!! "); - return 0; - } /* if */ - for (i = 0; i < m * n; i++) - flip[i] = matrix[i] * skal; - /* for */ - return flippy; - /* flippy = Free(flippy); */ -} /* extern double *JWDMMultMatSkalar */ - -/************************************************************************** - ROCKFLOW - Function: GetPriMatFromTriMat - - Task: - Gets a 9x9 Matrix from 3x3 matrix according to: - - a b c a b c - a b c d e f d e f - d e f ==> g h i g h i - g h i a b c a b c - d e f d e f - g h i g h i - - Parameter: (I: Input; R: Return; X: Both) - I: double *mat - - Return: - *mat2 - - Programming: - 07/2003 MB First Version -**************************************************************************/ - -int GetPriMatFromTriMat(double* mat1, double* mat_2) -{ - mat_2[0] = mat_2[3] = mat1[0]; - mat_2[1] = mat_2[4] = mat1[1]; - mat_2[2] = mat_2[5] = mat1[2]; - mat_2[6] = mat_2[9] = mat1[3]; - mat_2[7] = mat_2[10] = mat1[4]; - mat_2[8] = mat_2[11] = mat1[5]; - mat_2[12] = mat_2[15] = mat1[6]; - mat_2[13] = mat_2[16] = mat1[7]; - mat_2[14] = mat_2[17] = mat1[8]; - - mat_2[18] = mat_2[21] = mat1[0]; - mat_2[19] = mat_2[22] = mat1[1]; - mat_2[20] = mat_2[23] = mat1[2]; - mat_2[24] = mat_2[27] = mat1[3]; - mat_2[25] = mat_2[28] = mat1[4]; - mat_2[26] = mat_2[29] = mat1[5]; - mat_2[30] = mat_2[33] = mat1[6]; - mat_2[31] = mat_2[34] = mat1[7]; - mat_2[32] = mat_2[35] = mat1[8]; - - return 1; -} - -/************************************************************************** - ROCKFLOW - Funktion: MMultMatMat2 - Aufgabe: - Multiplikation zweier Matrizen der gleichen Grösse nach: - - a b x e f ==> axe bxf - c d g h cxg dxh - - Formalparameter: - E: *mat1, *mat2, m1, n1 - Ergebnis: - *matrix - -**************************************************************************/ -int MMultMatMat2(double* mat1, long m1, long n1, double* mat2, double* ergebnis) -{ - int i; - - for (i = 0; i < m1 * n1; i++) - ergebnis[i] = mat1[i] * mat2[i]; - - return 1; -} - -/************************************************************************** - ROCKFLOW - Function: TensorDrehDich - - Task: Dreht Tensor von von stromlinienorietierten in lokale physikalische - (Element) Koordinaten. - - Parameter: (I: Input; R: Return; X: Both) - I: double* d Tensor - double* velo Geschwindigkeitsvektor - - Return: - double* d - - Programming: - 08/2003 MB First Version, herausgelöst aus CalcEle3D - -**************************************************************************/ -double* TensorDrehDich(double* d, double* velo) - -{ - double zwa[18]; - double zwi[18]; - long l, ii; - double vg; - double fkt; - double zwo[9]; - - vg = MBtrgVec(velo, 3); - - if (d[0] > MKleinsteZahl || d[4] > MKleinsteZahl || d[8] > MKleinsteZahl) - /* Drehen: Stromrichtung - r,s,t */ - if (vg > MKleinsteZahl) - { - /* 1. Zeile */ - for (l = 0; l < 3; l++) - zwa[l] = velo[l]; - MNormiere(zwa, 3); - /* 2. Zeile */ - fkt = fabs(velo[0]); - ii = 0; - for (l = 1; l < 3; l++) - if (fabs(velo[l]) < fkt) - { - fkt = fabs(velo[l]); - ii = l; - } - zwo[ii] = 0.0; - zwo[(ii + 1) % 3] = velo[(ii + 2) % 3]; - zwo[(ii + 2) % 3] = -velo[(ii + 1) % 3]; - MNormiere(zwo, 3); - /* 3. Zeile */ - M3KreuzProdukt(zwa, zwo, zwi); - MNormiere(zwi, 3); - for (l = 0; l < 3; l++) - { - zwa[3 + l] = zwo[l]; - zwa[6 + l] = zwi[l]; - } - /* dreh^T * D * dreh */ - MMultMatMat(d, 3, 3, zwa, 3, 3, zwi, 3, 3); - MTranspoMat(zwa, 3, 3, zwo); - MMultMatMat(zwo, 3, 3, zwi, 3, 3, d, 3, 3); - } /* end if (vg > MKleinsteZahl) */ - /* end if (d[0] > MKleinsteZahl || d[4] > MKleinsteZahl || d[8] > MKleinsteZahl) */ - - return d; -} -#endif //#ifdef obsolete //05.03.2010 WW - /////////////////////////////////////////////////////////////////// // // Moved here by WW. 05.03.2010 diff --git a/FEM/mathlib.h b/FEM/mathlib.h index 94579a861..cc961e319 100644 --- a/FEM/mathlib.h +++ b/FEM/mathlib.h @@ -46,130 +46,6 @@ Mathematische Funktionen ########################################################################*/ -#ifdef obsolete // 01.2011 WW -extern int MGleichDouble(double zahl1, double zahl2, double tol); -/* MGleichDouble - Vergleicht zwei double-Zahlen unter - Beruecksichtigung einer Fehlertoleranz */ -extern int MOmega1D(double* vf, double r); -/* Berechnet 1D Ansatzfunktionen */ -extern int MOmega2D(double* vf, double r, double s); -/* Berechnet 2D Ansatzfunktionen */ -extern int MOmega3D(double* vf, double r, double s, double t); -/* Berechnet 3D Ansatzfunktionen */ -//#ifdef obsolete //WW. 11.2008 -extern int MOmega3DTetrahedron(double* vf, double xx, double yy, double zz, long number); -extern int MOmega2DTriangle(double* vf, double r, double s, long number); -//#endif -/* Berechnet 2D Ansatzfunktionen fuer Dreiecke */ -// extern int MPhi2D(double *vf,double r, double s); -/* Berechnet 2D Testfunktionen */ -extern int MPhi2D_SUPG(double* vf, double r, double s, double* alpha); -/* Berechnet 2D Testfunktionen (SUPG) */ -extern int MPhi3D(double* vf, double r, double s, double t); -/* Berechnet 3D Testfunktionen */ -extern int MPhi3D_SUPG(double* vf, double r, double s, double t, double* alpha); -/* Berechnet 3D Testfunktionen (SUPG)*/ -extern int MGradOmega2D(double* vf, double r, double s); -/* Berechnet Gradient der 2D Ansatzfunktionen */ -extern int MGradOmega3D(double* vf, double r, double s, double t); -/* Berechnet Gradient der 3D Ansatzfunktionen */ -extern int MGradPhi2D(double* vf, double r, double s); -/* Berechnet Gradient der 2D Testfunktionen */ -extern int MGradPhi3D(double* vf, double r, double s, double t); -/* Faktoren fuer die X Punkt Gauss-Integration */ -extern void MGetCoor(int typ, long j, double* r, double* s, double* t); -/* lokale Koordinaten der Eckpunkte */ - -/* MBtrgVec - Betrag von Vektor */ - -/* MAngleVectors - Winkel zwischen Vektoren */ /* MB */ -double MAngleVectors(double* v1, double* v2); -/* MNormiere - Normiert Vektoren */ -void MNormiere(double* vec, long n); -/* M2Determinante - Determinante einer 2x2Matrix */ -extern double M2Determinante(double* matrix); -/* M3Determinante - Determinante einer 3x3Matrix */ -// extern double M3Determinante ( double *m ); -/* Mxg2Determinante - Determinante einer beliebigen Matrix - mit goesseren Ausmassen als 2x2 */ -/* M4Determinante - Determinante einer 4x4Matrix */ -// extern double M4Determinante ( double *m ); -extern double Mxg2Determinante(double* matrix, long m, long n); -/* MTranspoVec - Transponieren beliebiger Vektoren */ -extern void MTranspoVec(double* vec, long g); -/* MTranspoMat - Transponieren beliebiger Matrizen */ -extern void MTranspoMat(double* mat1, long m, long n, double* mat2); -/* M2Invertiere - Invertiert 2x2 Matrizen */ -extern void M2Invertiere(double* m); -/* M2InvertiereUndTransponiere - Invertiert und transponiert 2x2 Matrizen */ -extern void M2InvertiereUndTransponiere(double* m); -/* M3Invertiere - Invertiert 3x3 Matrizen */ -extern void M3Invertiere(double* m); -/* MInvertiere - Invertieren beliebier regulaerer Matrizen */ -extern void MInvertiere(double* matrix, long m, long n); -/* MAddVektoren - Addition zweier beliebier Vektoren */ -extern int MAddVektoren(double* v1, double* v2, double* vout, long g); -/* MAddSkalVektoren - Vektoren mit Skalar multiplizieren und dann addieren */ -// WW extern int MAddSkalVektoren(double *v1, double m1, double *v2, double m2, double *vout, long g); -/* MAddMatrizen - Addition zweier beliebier Matrizen */ -extern int MAddMatrizen(double* m1, double* m2, double* mout, long m, long n); -/* MMultVecSkalar - Multiplikation Vektor mit Skalarwert */ -extern int MMultVecSkalar(double* vec, double skal, long g); -/* MMultMatSkalar - Multiplikation Matrix mit Skalarwert */ -extern int MMultMatSkalar(double* matrix, double skal, long m, long n); -/* MSkalarprodukt - Skalarprodukt zweier beliebiger Vektoren */ -// extern double MSkalarprodukt ( double *vec1, double *vec2, long g ); -/* M3KreuzProdukt - Kreuzprodukt von 3D-Vektoren */ -// extern int M3KreuzProdukt( double *vec1, double *vec2 , double *vec); -/*########################################################################## - Auf Matrizen und Vektoren aufbauende Funktionen - ########################################################################*/ - -int MKTF2Dr2D(double* vec1, double winkel, double* vec2); -int MKTFMat2Dr2D(double winkel, double* tmat); -int MKTF2Dt2D(double* vec1, double dx, double dy, double* vec2); -int MKTF3Dr2D(double* vec1, double* vec2, double winkel, double* vec); -void MKTFMat3Dr2D(double* vec1, double* vec2, double winkel, double* mat); -/* was ist das alles ??? (msr) */ - -/*########################################################################## - Prueffunktion fuer Matrix - ########################################################################*/ -/* MBistDuDiagMat - Prueft auf Diagonalitaet einer Matrix */ -extern int MBistDuDiagMat(double* matrix, long m, long n); - -/*########################################################################## - Sortierfunktionen - ########################################################################*/ -extern void MQSort_LongDouble(void* DataSets, const int NumberOfDataSets, const int SiezeOfDataSet); -extern int MCompare_for_MQSort_LongDouble(const void* arg1, const void* arg2); - -extern double* TensorDrehDich(double* d, double* velo); -extern int GetPriMatFromTriMat(double* mat1, double* mat2); -extern int MMultMatMat2(double* mat1, long m1, long n1, double* mat2, double* ergebnis); - -/*########################################################################## - eventuell noch nuetzliche Funktionen, - die im Moment nicht gebraucht werden - ########################################################################*/ -extern double* JWDMMultMatSkalar(double* matrix, double skal, long m, long n); - -// extern long IsNodeInsideTriangle (long n1, long n2, long n3, long node); -// CC -extern double CalcTriangleArea(long n1, long n2, long n3); -// extern long IsNodeInsidePlain (long n1, long n2, long n3, long node); - -extern int MOmega2D_9N(double* vf, double r, double s); -extern int MPhi2D_9N(double* vf, double r, double s); -extern int MGradOmega2D_9N(double* vf, double r, double s); -extern int MGradPhi2D_9N(double* vf, double r, double s); - -extern int MOmega3D_20N(double* vf, double r, double s, double t); -extern int MPhi3D_20N(double* vf, double r, double s, double t); -extern int MGradOmega3D_20N(double* vf, double r, double s, double t); -extern int MGradPhi3D_20N(double* vf, double r, double s, double t); -#endif // 05.03.2010. WW //#ifdef obsolete - /* Ermittelt min */ extern double MMin(double, double); /* Ermittelt max */