src/expde/formulas/boundy.cc

Go to the documentation of this file.

//    expde: expression templates for partial differential equations.
//    Copyright (C) 2001  Christoph Pflaum
//    This program is free software; you can redistribute it and/or modify
//    it under the terms of the GNU General Public License as published by
//    the Free Software Foundation; either version 2 of the License, or
//    (at your option) any later version.
//
//    This program is distributed in the hope that it will be useful,
//    but WITHOUT ANY WARRANTY; without even the implied warranty of
//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//    GNU General Public License for more details.
//
//                 SEE  Notice1.doc made by 
//                 LAWRENCE LIVERMORE NATIONAL LABORATORY
//

// ------------------------------------------------------------
//
// boundary.cc
//
// ------------------------------------------------------------

// Include level 0:
#ifdef COMP_GNUOLD
 #include <iostream.h>
 #include <fstream.h>
 #include <math.h>
#else
 #include <iostream>
 #include <fstream>
 #include <cmath>
#endif 

// Include level 1:
#include "../paramete.h"
#include "../abbrevi.h"
#include "../math_lib/math_lib.h"

// Include level 4:
#include "../formulas/boundy.h"
#include "subdiv.h"




// now in math_lib.h

#define x0 V0.x
#define y0 V0.y
#define z0 V0.z

#define x1 V1.x
#define y1 V1.y
#define z1 V1.z

#define x2 V2.x
#define y2 V2.y
#define z2 V2.z

#define x3 V3.x
#define y3 V3.y
#define z3 V3.z


/*
inline double Calc_det(const D3vector& V0, const D3vector& V1,
                       const D3vector& V2, const D3vector& V3) {
  return (x1*y2*z0 - 
          x1*y3*z0 - x0*y2*z1 + x0*y3*z1 - x1*y0*z2 + x0*y1*z2 - x0*y3*z2 + 
          x1*y3*z2 + x3*(y1*z0 - y2*z0 - y0*z1 + y2*z1 + y0*z2 - y1*z2) + 
          x1*y0*z3 - x0*y1*z3 + x0*y2*z3 - x1*y2*z3 + 
          x2*(y3*z0 + y0*z1 - y3*z1 - y0*z3 + y1*(-z0 + z3)));
}
*/

#define Wx (V0.x - V1.x)
#define Wy (V0.y - V1.y)
#define Wz (V0.z - V1.z)

#define Vx (V0.x - V2.x)
#define Vy (V0.y - V2.y)
#define Vz (V0.z - V2.z)

#define Vcrossx (Wy * Vz - Wz * Vy)
#define Vcrossy (Wx * Vz - Wz * Vx)
#define Vcrossz (Wy * Vx - Wx * Vy)

inline double Calc_det_triangle(const D3vector& V0, const D3vector& V1,
                                const D3vector& V2) {
  return sqrt(Vcrossx*Vcrossx + Vcrossy*Vcrossy + Vcrossz*Vcrossz);
}

void BoCeData::Set_h(double *h, Bo_description& desc) {
  int num;

  for(num=0;num<number_points;++num) {
    if(edge_point(num)==edge_poi_typ) {
      h[num] = desc.h(corner(num),edge_dir(num));
    }
  }


}

void BoCeData::Calc_FE_elements(Bo_description& desc) {
  Calculator_FE calculator;
  Tetraeder_storage* tet;
  double h[24];
  double deter;
  D3vector V0;
  D3vector V1;
  D3vector V2;
  D3vector V3;
  int num;

  // 1. Initialisierung: Informationen ueber Randzelle abspeichern
  Initialize(desc);

  // 2. Berechnung der Tetraeder
  calculator.Start(&desc,this);
  Act_on_subdivision(&desc,&calculator);

  // An letzten inneren Tet: Tets am Rand haengen.
  if(Tets!=NULL) {
    tet=Tets;
    if(tet->Next()!=NULL) {
      while(tet->Next()!=NULL) { tet=tet->Next(); }
    }
    tet->Put_Next(Tets_boundary);
  }
  else Tets = Tets_boundary;


  // 3. Berechnung der Determinante und lokaler Matrix der Tetraeder
  // Berechnung der Maschenweiten
  hhh = new double[number_points];
  for(num=0;num<number_points;++num) {
    if(edge_point(num)==edge_poi_typ) {
      h[num] = desc.h(corner(num),edge_dir(num));
      hhh[num] = h[num];
    }
  }


  


  
  for(tet=Tets;tet!=NULL;tet=tet->Next()) {
    // Berechnung der Koordinaten
    V0 = coord(tet->N0(),h);
    V1 = coord(tet->N1(),h);
    V2 = coord(tet->N2(),h);
    V3 = coord(tet->N3(),h);

    // Determinante
    deter = Calc_det(V0,V1,V2,V3);
    tet->Put_det(deter);
    if(developer_version) {
      if(deter < 0.0)
        cout << "\n Negatives Gewicht in Tet! " << endl;
    }
    // Speicher fuer lokale Matrix
    tet->Xm0 = y3 * (z1 - z2) + y1 * (z2 - z3) + y2 * (z3 - z1);
    tet->Xm1 = y3 * (z2 - z0) + y2 * (z0 - z3) + y0 * (z3 - z2);
    tet->Xm2 = y3 * (z0 - z1) + y0 * (z1 - z3) + y1 * (z3 - z0);
    tet->Xm3 = y2 * (z1 - z0) + y1 * (z0 - z2) + y0 * (z2 - z1);
    
    tet->Ym0 = x3 * (z2 - z1) + x2 * (z1 - z3) + x1 * (z3 - z2);
    tet->Ym1 = x3 * (z0 - z2) + x0 * (z2 - z3) + x2 * (z3 - z0);
    tet->Ym2 = x3 * (z1 - z0) + x1 * (z0 - z3) + x0 * (z3 - z1);
    tet->Ym3 = x2 * (z0 - z1) + x0 * (z1 - z2) + x1 * (z2 - z0);

    tet->Zm0 = x3 * (y1 - y2) + x1 * (y2 - y3) + x2 * (y3 - y1);
    tet->Zm1 = x3 * (y2 - y0) + x2 * (y0 - y3) + x0 * (y3 - y2);
    tet->Zm2 = x3 * (y0 - y1) + x0 * (y1 - y3) + x1 * (y3 - y0);
    tet->Zm3 = x2 * (y1 - y0) + x1 * (y0 - y2) + x0 * (y2 - y1);

    // Determinante ueber Randflaechen
    if(tet->Boundary_tet()) {
      tet->Put_det_surface(Calc_det_triangle(V0,V1,V2));
    }
  }
  // 4. Speicher fuer Liste der Zeiger auf Variablen initialisieren 
  if(Exists_bocellpoint())
    vars   = new double*[number_points+1];
  else
    vars   = new double*[number_points];

  // 5. Only for looking for a mistake
  // Analyse_angles(desc);
}

void BoCeData::Analyse_angles(Bo_description& desc) {
  Tetraeder_storage* tet;
  bool mistake;
  mistake = false;
  for(tet=Tets;tet!=NULL;tet=tet->Next()) {
    if(tet->Check_angles(&desc)==false) mistake=true;
  }
  if(mistake) {
    Print();
    for(tet=Tets;tet!=NULL;tet=tet->Next()) {
      if(tet->Check_angles(&desc)==false) {
        cout << "\n Corners of Tet with bad angles: " <<
          tet->N0() << ", " << 
          tet->N1() << ", " << 
          tet->N2() << ", " << 
          tet->N3() << ", " << endl;
      }
    }
    cout << "Meshsize: " << Meshsize() << endl;
  }
}

// Funktionen zum Testen:  muss spaeter wohl weg:
void Bo_description::Test_Initialisierung() {
  H = 1.0;

 // Beispiel fuer einen Tetraeder bei der Ecke: WSDd
  /*
 cout << " Beispiel fuer einen Tetraeder bei der Ecke: WSDd \n";
 Put_corner_point(WSDd);
 Put_edge_point(WSDd,Edir,0.5);
 Put_edge_point(WSDd,Ndir,0.5);
 Put_edge_point(WSDd,Tdir,0.5);
  */
  
  /*
  // Beispiel fuer eine Kante im Gebiet
  cout << " Beispiel fuer eine Kante im Gebiet  \n";
  Put_corner_point(WSDd);
  Put_corner_point(WSTd);

  Put_edge_point(WSDd,Edir,0.5);
  Put_edge_point(WSDd,Ndir,0.5);

  Put_edge_point(WSTd,Edir,0.5);
  Put_edge_point(WSTd,Ndir,0.5);
  
  Put_edge(SWed);
  */

  // Beispiel fuer einen halben Wuerfel
  cout << " Beispiel fuer einen halben Wuerfel  \n";
  Put_corner_point(WSDd);
  Put_corner_point(WNDd);
  Put_corner_point(ESDd);
  Put_corner_point(ENDd);

  Put_edge_point(WSDd,Tdir,0.5);
  Put_edge_point(WNDd,Tdir,0.5);
  Put_edge_point(ESDd,Tdir,0.5);
  Put_edge_point(ENDd,Tdir,0.5);
  
  Put_edge(SDed);
  Put_edge(NDed);
  Put_edge(EDed);
  Put_edge(WDed);


  // Beispiel fuer zwei nicht zusammenhaengende Tetraeder
  /*
  cout << " Beispiel fuer zwei nicht zusammenhaengende Tetraeder \n";
  Put_corner_point(WSDd);
  Put_edge_point(WSDd,Edir,0.5);
  Put_edge_point(WSDd,Ndir,0.5);
  Put_edge_point(WSDd,Tdir,0.5);

  Put_corner_point(WNDd);
  Put_edge_point(WNDd,Edir,0.5);
  Put_edge_point(WNDd,Sdir,0.5);
  Put_edge_point(WNDd,Tdir,0.5);
*/
}


void BoCeData::Initialize(Bo_description& desc) {
  static int i,k,number;

  bocellpoint=NULL; 

  number=0;
  // 1. Anzahl der Eckpunkte und Kantenpunkte zaehlen
  for(i=0;i<8;++i)  if(desc.corner_point((dir_sons)i)) ++number;
  for(i=0;i<8;++i) for(k=0;k<3;++k) 
    if(desc.edge_point((dir_sons)i,(dir_3D)(k<<1))) ++number;
  // 2. Speicher bereitstellen
  Allocate(number);

  // 3. Informationen setzen
  desc.Set_num_cellpoi(number_points);
  number=0;
  for(i=0;i<8;++i)  if(desc.corner_point((dir_sons)i)) {
    Put_corner_point(number,(dir_sons)i);
    desc.Set_num_corner((dir_sons)i,number);
    ++number;
  }
  for(i=0;i<8;++i) {
    if(desc.edge_point((dir_sons)i,Wdir)) {
      Put_edge_point(number,(dir_sons)i,(dir_3D)((~i)&1));
      desc.Set_num_edge_poi((dir_sons)i,(dir_3D)((~i)&1),number);
      ++number;
    }
    if(desc.edge_point((dir_sons)i,Sdir)) {
      Put_edge_point(number,(dir_sons)i,(dir_3D)(2+(((~i)>>1)&1)));
      desc.Set_num_edge_poi((dir_sons)i,(dir_3D)(2+(((~i)>>1)&1)),number);
      ++number;
    }
    if(desc.edge_point((dir_sons)i,Ddir)) {
      Put_edge_point(number,(dir_sons)i,(dir_3D)(4+(((~i)>>2)&1)));
      desc.Set_num_edge_poi((dir_sons)i,(dir_3D)(4+(((~i)>>2)&1)),number);
      ++number;
    }
  }
  // 4. Maschenweite setzen
  Put_Meshsize(desc.Meshsize());
  // 5. Kanteninfo setzen
  Put_edge_info(desc.edge_info());

  // 5. Setz Zeiger auf Bocelldata
  desc.Put_bocedata(this);
}



// Funktion, die die Informationen der Randzelle ausdruckt
void BoCeData::Print() {
  int i;
  cout << "\n 1.Teil: Punkte:";
  for(i=0;i<number_points;++i) {
    cout << "\n Nummer " << i;
    if(edge_point(i)) {
      cout << " ist Kantenpunkt bei Ecke: ";
      Print(corner(i));
      cout << " in Richtung: ";
      Print(edge_dir(i));
     }
    else {
      cout << " ist Eckpunkt    bei Ecke: ";
      Print(corner(i));
      cout << "                  ";
    }
  }
  cout << "\n 2.Teil: Kanten: \n ";
  for(i=0;i<12;++i) if(Edge((Edges_cell)i)) {
    cout << "Kante im Gebiet: ";
    Print((Edges_cell)i);
    cout << endl;
  }
  cout << endl;

  cout << " Determinante of tets: " << endl; 
  Tetraeder_storage* tet;
  for(tet=Give_tets();tet!=NULL;tet=tet->Next()) {
    cout << " Tet: det: "  << tet->Det() << endl;
  }
  cout << " Connection of tets: " << endl; 
  for(tet=Give_tets();tet!=NULL;tet=tet->Next()) {
    cout << " new tet num: " << endl;
    cout << tet->N0() << ", " << tet->N1() << ", " 
         << tet->N2() << ", " << tet->N3() << endl;
  }
}

// Ganz einfache Hilfs-Funktionen  
void BoCeData::Print(dir_3D d) {
  if(d==Wdir) cout << "Wdir"; 
  if(d==Edir) cout << "Edir"; 
  if(d==Sdir) cout << "Sdir"; 
  if(d==Ndir) cout << "Ndir"; 
  if(d==Ddir) cout << "Ddir"; 
  if(d==Tdir) cout << "Tdir"; 
}


void BoCeData::Print(dir_sons d) {
  if(d==WSDd) cout << "WSDd"; 
  if(d==ESDd) cout << "ESDd"; 
  if(d==WNDd) cout << "WNDd"; 
  if(d==ENDd) cout << "ENDd"; 
  if(d==WSTd) cout << "WSTd"; 
  if(d==ESTd) cout << "ESTd"; 
  if(d==WNTd) cout << "WNTd"; 
  if(d==ENTd) cout << "ENTd"; 
}


void BoCeData::Print(Edges_cell ed) {
  if(ed==SDed) cout << "SDed";
  if(ed==NDed) cout << "NDed";
  if(ed==STed) cout << "STed";
  if(ed==NTed) cout << "NTed";
  if(ed==WDed) cout << "WDed";
  if(ed==EDed) cout << "EDed";
  if(ed==WTed) cout << "WTed";
  if(ed==ETed) cout << "ETed";
  if(ed==SWed) cout << "SWed";
  if(ed==SEed) cout << "SEed";
  if(ed==NWed) cout << "NWed";
  if(ed==NEed) cout << "NEed";
}

---