src/expde/formulas/subdiv.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
//

// ------------------------------------------------------------
//
// subdiv.cc
// ruft spaeter subdivI.h subdivII.h subdivIII.h subdivIV.h auf
//
// ------------------------------------------------------------

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


// Include level 0:
#include "../parser.h"

// 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"



// Funktion die eine Randzelle in 
// Tetraeder zerlegt

/* 
    Function which calculates angle angle(P1,P3,P2) 
    P1 x            x  P2
       |            |
    s1 |            |  s2
       |            |
       x-----x------x
      P4     P3     P5
          s3    s4

*/

double calc_special_angle(double s1, double  s3, double s2, double s4) {
  return 180.0 - (my_atan(s1/s3) + my_atan(s2/s4)) / M_PI * 180.0; 
}

// devisions for prism:
// these functions give back the type of the face
// (P_5,P_2,P_1,P_4) at the point P_1
                   
// Zerlegung eines Prismas in Tetraeder
// Nimm immer 1. Version
f_type Act_on_Prism(Bo_description* Description, Calculator_FE* Act,
                  f_type x_face_type, f_type y_face_type,
                  Edge_Corner_point poi_0, Edge_Corner_point poi_1, 
                  Edge_Corner_point poi_2, Edge_Corner_point poi_3,
                  Edge_Corner_point poi_4, Edge_Corner_point poi_5) {
  if(x_face_type==an && y_face_type==an) {
    Act->Tetraeder(poi_0,poi_2,poi_3,poi_1);
    Act->Tetraeder(poi_2,poi_3,poi_1,poi_4);
    Act->Tetraeder(poi_5,poi_3,poi_2,poi_4);
    return an;
  }
  if(x_face_type==spitz && y_face_type==an) {
    Act->Tetraeder(poi_2,poi_0,poi_1,poi_4);
    Act->Tetraeder(poi_3,poi_0,poi_2,poi_4);
    Act->Tetraeder(poi_5,poi_3,poi_2,poi_4);
    return an;
  }
  if(x_face_type==an && y_face_type==spitz) {
    Act->Tetraeder(poi_5,poi_3,poi_0,poi_1);
    Act->Tetraeder(poi_2,poi_5,poi_0,poi_1);
    Act->Tetraeder(poi_4,poi_5,poi_1,poi_3);
    return spitz;
  }
  //  if(x_face_type==spitz && y_face_type==spitz) {
  Act->Tetraeder(poi_5,poi_2,poi_1,poi_0);
  Act->Tetraeder(poi_5,poi_3,poi_0,poi_4);
  Act->Tetraeder(poi_1,poi_5,poi_0,poi_4);    
  return spitz;
  //  }
}

// Zerlegung eines Prismas in Tetraeder
// Nimm immer beste Version
f_type Act_on_Prism_best(Bo_description* Description, Calculator_FE* Act,
                  f_type x_face_type, f_type y_face_type,
                  Edge_Corner_point poi_0, Edge_Corner_point poi_1, 
                  Edge_Corner_point poi_2, Edge_Corner_point poi_3,
                  Edge_Corner_point poi_4, Edge_Corner_point poi_5) {
  if(developer_version) {
    if(poi_2.edge_point!=edge_poi_typ ||
       poi_1.edge_point!=edge_poi_typ ||
       poi_4.edge_point!=edge_poi_typ ||
       poi_5.edge_point!=edge_poi_typ)
      cout << " \n Mistake in Act_on_Prism_best ! \n ";
  }

  if(x_face_type==an && y_face_type==an) {
    if(parallel_version)
      Tell_about_round_of_problem(Description->h(poi_2.corner, poi_2.d),
                                  Description->h(poi_1.corner, poi_1.d));
    if(Description->h(poi_2.corner, poi_2.d) >
       Description->h(poi_1.corner, poi_1.d)) {
      Act->Tetraeder(poi_0,poi_2,poi_3,poi_1);
      Act->Tetraeder(poi_2,poi_3,poi_1,poi_4);
      Act->Tetraeder(poi_5,poi_3,poi_2,poi_4);
      return an;
    }
    else  {
      Act->Tetraeder(poi_0,poi_2,poi_3,poi_1);
      Act->Tetraeder(poi_2,poi_3,poi_1,poi_5);
      Act->Tetraeder(poi_5,poi_3,poi_1,poi_4);
      return spitz;
    }
  }
  if(x_face_type==spitz && y_face_type==an) {
    Act->Tetraeder(poi_2,poi_0,poi_1,poi_4);
    Act->Tetraeder(poi_3,poi_0,poi_2,poi_4);
    Act->Tetraeder(poi_5,poi_3,poi_2,poi_4);
    return an;
  }
  if(x_face_type==an && y_face_type==spitz) {
    Act->Tetraeder(poi_5,poi_3,poi_0,poi_1);
    Act->Tetraeder(poi_2,poi_5,poi_0,poi_1);
    Act->Tetraeder(poi_4,poi_5,poi_1,poi_3);
    return spitz;
  }
  //  if(x_face_type==spitz && y_face_type==spitz) {
  if(parallel_version)
      Tell_about_round_of_problem(Description->h(poi_5.corner, poi_5.d),
                                  Description->h(poi_4.corner, poi_4.d));
  if(Description->h(poi_5.corner, poi_5.d) >
     Description->h(poi_4.corner, poi_4.d)) {
    Act->Tetraeder(poi_5,poi_2,poi_1,poi_0);
    Act->Tetraeder(poi_5,poi_3,poi_0,poi_4);
    Act->Tetraeder(poi_1,poi_5,poi_0,poi_4);    
    return spitz;
  }
  else {
    Act->Tetraeder(poi_5,poi_2,poi_4,poi_0);
    Act->Tetraeder(poi_5,poi_3,poi_0,poi_4);
    Act->Tetraeder(poi_1,poi_2,poi_0,poi_4);    
    return an;
  }


  //  }
}

// Zerlegung eines Prismas in Tetraeder
// Take version which is better for |01| < |02|
//                          and for |34| < |35|
f_type Act_on_Prism_short(Bo_description* Description, Calculator_FE* Act,
                  f_type x_face_type, f_type y_face_type,
                  Edge_Corner_point poi_0, Edge_Corner_point poi_1, 
                  Edge_Corner_point poi_2, Edge_Corner_point poi_3,
                  Edge_Corner_point poi_4, Edge_Corner_point poi_5) {
  if(x_face_type==an && y_face_type==an) {
    Act->Tetraeder(poi_0,poi_2,poi_3,poi_1);
    Act->Tetraeder(poi_2,poi_3,poi_1,poi_4);
    Act->Tetraeder(poi_5,poi_3,poi_2,poi_4);
    return an;
  }
  if(x_face_type==spitz && y_face_type==an) {
    Act->Tetraeder(poi_2,poi_0,poi_1,poi_4);
    Act->Tetraeder(poi_3,poi_0,poi_2,poi_4);
    Act->Tetraeder(poi_5,poi_3,poi_2,poi_4);
    return an;
  }
  if(x_face_type==an && y_face_type==spitz) {
    Act->Tetraeder(poi_5,poi_3,poi_0,poi_1);
    Act->Tetraeder(poi_2,poi_5,poi_0,poi_1);
    Act->Tetraeder(poi_4,poi_5,poi_1,poi_3);
    return spitz;
  }
  //  if(x_face_type==spitz && y_face_type==spitz) {
  Act->Tetraeder(poi_5,poi_2,poi_1,poi_0);
  Act->Tetraeder(poi_5,poi_3,poi_0,poi_4);
  Act->Tetraeder(poi_1,poi_5,poi_0,poi_4);    
  return spitz;
  //  }
}

// Zerlegung eines Prismas in Tetraeder
// Take version which is better for |01| > |02|
//                          and for |34| > |35|
f_type Act_on_Prism_long(Bo_description* Description, Calculator_FE* Act,
                  f_type x_face_type, f_type y_face_type,
                  Edge_Corner_point poi_0, Edge_Corner_point poi_1, 
                  Edge_Corner_point poi_2, Edge_Corner_point poi_3,
                  Edge_Corner_point poi_4, Edge_Corner_point poi_5) {
  if(x_face_type==an && y_face_type==an) {
    Act->Tetraeder(poi_0,poi_2,poi_3,poi_1);
    Act->Tetraeder(poi_2,poi_3,poi_1,poi_5);
    Act->Tetraeder(poi_5,poi_3,poi_1,poi_4);
    return spitz;
  }
  if(x_face_type==spitz && y_face_type==an) {
    Act->Tetraeder(poi_2,poi_0,poi_1,poi_4);
    Act->Tetraeder(poi_3,poi_0,poi_2,poi_4);
    Act->Tetraeder(poi_5,poi_3,poi_2,poi_4);
    return an;
  }
  if(x_face_type==an && y_face_type==spitz) {
    Act->Tetraeder(poi_5,poi_3,poi_0,poi_1);
    Act->Tetraeder(poi_2,poi_5,poi_0,poi_1);
    Act->Tetraeder(poi_4,poi_5,poi_1,poi_3);
    return spitz;
  }
  //  if(x_face_type==spitz && y_face_type==spitz) {
  Act->Tetraeder(poi_5,poi_2,poi_4,poi_0);
  Act->Tetraeder(poi_5,poi_3,poi_0,poi_4);
  Act->Tetraeder(poi_1,poi_2,poi_0,poi_4);    
  return an;
  //  }
}


#include "subdivI.h"
#include "subdivII.h"
#include "subdiIII.h"
#include "subdivIV.h"

// Zerlegung der Randzelle
void Act_on_subdivision(Bo_description* Description, Calculator_FE* Act) {
  int i;
  bool corners[8];
  bool stop;

  // Vorbereitung Infos ueber Ecken setzen:
  // Wichtig fuer doppelte Randzellen!
  for(i=0;i<8;++i) 
    corners[i] = Description->corner_point((dir_sons)i);

  Act_on_subdivision_partI(Description,Act,corners);
  // Entscheidung ob noch weitere Ecken behandelt werden muessen
  stop = true;
  for(i=0;i<8;++i) if(corners[i]) stop = false;
  Act_on_subdivision_partII(Description,Act,corners,stop);
  Act_on_subdivision_partIII(Description,Act,corners,stop);
  Act_on_subdivision_partIV(Description,Act,corners,stop);
};


//  Memberfunktion


void Calculator_FE::Tetraeder(Edge_Corner_point p0, Edge_Corner_point p1,
                              Edge_Corner_point p2, Edge_Corner_point p3) { 
  if(p0.edge_point!=edge_poi_typ &&
     p1.edge_point==edge_poi_typ &&
     p2.edge_point==edge_poi_typ &&
     p3.edge_point==edge_poi_typ) {
    bo_cell->Add_Tetraeder_bo(bo_desc,
                              bo_desc->Num(p2),bo_desc->Num(p1),
                              bo_desc->Num(p3),bo_desc->Num(p0));
  }
  else if(p0.edge_point==edge_poi_typ &&
          p1.edge_point!=edge_poi_typ &&
          p2.edge_point==edge_poi_typ &&
          p3.edge_point==edge_poi_typ) {
    bo_cell->Add_Tetraeder_bo(bo_desc,
                              bo_desc->Num(p2),bo_desc->Num(p3),
                              bo_desc->Num(p0),bo_desc->Num(p1));
  }
  else if(p0.edge_point==edge_poi_typ &&
          p1.edge_point==edge_poi_typ &&
          p2.edge_point!=edge_poi_typ &&
          p3.edge_point==edge_poi_typ) {
    bo_cell->Add_Tetraeder_bo(bo_desc,
                              bo_desc->Num(p0),bo_desc->Num(p3),
                              bo_desc->Num(p1),bo_desc->Num(p2));
  }
  else if(p0.edge_point==edge_poi_typ &&
          p1.edge_point==edge_poi_typ &&
          p2.edge_point==edge_poi_typ &&
          p3.edge_point!=edge_poi_typ) {
    bo_cell->Add_Tetraeder_bo(bo_desc,
                              bo_desc->Num(p0),bo_desc->Num(p1),
                              bo_desc->Num(p2),bo_desc->Num(p3));
  }
  else
    bo_cell->Add_Tetraeder(bo_desc,
                           bo_desc->Num(p0),bo_desc->Num(p1),
                           bo_desc->Num(p2),bo_desc->Num(p3));
}

---