src/expde/formulas/subdivI.h

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

// ------------------------------------------------------------
//
// subdivI.h
//
// ------------------------------------------------------------

#ifndef SUBDIVI_H_
#define SUBDIVI_H_

void Act_on_subdivision_partI(Bo_description* Description, 
                              Calculator_FE* Act, bool* corners) {
  int i, k;

  Edge_Corner_point poi_0,   poi_1,   poi_2,   poi_3,   poi_4;
  Edge_Corner_point poi_5,   poi_6,   poi_7,   poi_8;
  Edge_Corner_point poi_A,   poi_B,   poi_C,   poi_D;
  Edge_Corner_point poi_A2,  poi_B2;
  Edge_Corner_point poi_boA, poi_boB;
  Edge_Corner_point poi_boA2, poi_boB2;
  Edge_Corner_point poi_oppA,   poi_oppB;

  dir_sons    oppA, oppB, corner;
  dir_3D    dirA, dirB, dirC;
  bool stop;
  int ori;

  f_type x_face_type;
  f_type y_face_type;

  f_type A_face_type, B_face_type;

  // leere Initialisierung
  ori = 0;

  stop = false;
  // stop = true;
  //  Boundary cell: 3.2.
  // ---------------------
  //  Ganz einfache Tetraeder in der Ecke
  poi_4.edge_point = corner_poi_typ;
  poi_1.edge_point = edge_poi_typ;
  poi_2.edge_point = edge_poi_typ;
  poi_3.edge_point = edge_poi_typ;

  for(i=0;i<8 && !stop;++i) {
    poi_1.corner = poi_2.corner = poi_3.corner = poi_4.corner = (dir_sons)i;
    poi_1.d = opposite3D(xCoord((dir_sons)i));
    poi_2.d = opposite3D(yCoord((dir_sons)i));
    poi_3.d = opposite3D(zCoord((dir_sons)i));

    if(Description->edge_point(poi_1) &&
       Description->edge_point(poi_2) &&
       Description->edge_point(poi_3))    {
      corners[i] = false;
      // stop=true;
      if((((i&1) + ((i>>1)&1) + ((i>>2)&1))&1)==0) { // Orientierung 
        //  1. Version 
        Act->Tetraeder(poi_1,poi_4,poi_2,poi_3);
      }
      else {
        // 2. Version
        Act->Tetraeder(poi_1,poi_4,poi_3,poi_2);
      }
    }
  }

  //  stop = false;
  //  stop = true;
  //  Boundary cell : 3.1.
  // ----------------------
  // III. Ein Prisma
  for(i=0;i<12 && !stop;++i) {
    // Ganze Kante im Gebiet
    if(Description->Edge((Edges_cell)i)) {
      // Zwei Ecken der Kante
      poi_A.edge_point = poi_B.edge_point = corner_poi_typ;

      poi_A.corner = Transform((Edges_cell)i,Ld);
      poi_B.corner = Transform((Edges_cell)i,Rd);

      // Richtungen
      dirC = Trans_dir((Edges_cell)i,Rd);
      dirA = OrthoA((Edges_cell)i);
      dirB = OrthoB((Edges_cell)i);

      // Vier Kantepunkte
      poi_boA.edge_point  = poi_boB.edge_point  = edge_poi_typ;
      poi_boA2.edge_point = poi_boB2.edge_point = edge_poi_typ;

      // dL: A
      poi_boA.corner = poi_A.corner;
      poi_boA.d = opposite3D(OrthoB((Edges_cell)i));
      poi_boA2.corner = poi_A.corner;
      poi_boA2.d = opposite3D(OrthoA((Edges_cell)i));
      // dR: B
      poi_boB.corner = poi_B.corner;
      poi_boB.d = opposite3D(OrthoB((Edges_cell)i));
      poi_boB2.corner = poi_B.corner;
      poi_boB2.d = opposite3D(OrthoA((Edges_cell)i));

      if(Description->edge_point(poi_boA) &&
         Description->edge_point(poi_boA2) &&
         Description->edge_point(poi_boB) &&
         Description->edge_point(poi_boB2)) {
        corners[poi_A.corner] = false;
        corners[poi_B.corner] = false;

        //      stop = true;
        // Orientierung:
        if(dirC==Wdir) ori=0;
        if(dirC==Edir) ori=0;
        if(dirC==Ndir) ori=1;
        if(dirC==Sdir) ori=1;
        if(dirC==Tdir) ori=0;
        if(dirC==Ddir) ori=0;
        
        
        if(((OneDpart(dirA)+OneDpart(dirB))&1)==ori) {
          poi_1 = poi_boA;        poi_2 = poi_boA2;
          poi_0 = poi_A;          poi_4 = poi_boB;
          poi_5 = poi_boB2;       poi_3 = poi_B;
        }
        else {
          poi_1 = poi_boA2;       poi_2 = poi_boA;
          poi_0 = poi_A;          poi_4 = poi_boB2;
          poi_5 = poi_boB;        poi_3 = poi_B;

        }
        dirC = direction_from_to(poi_0.corner,poi_3.corner);
        dirA = poi_1.d;
        dirB = poi_2.d;
        x_face_type = Calc_f_type(dirA,dirC);
        y_face_type = Calc_f_type(dirB,dirC);
        //      Act_on_Prism(Description,Act,x_face_type,y_face_type,
        //                   poi_0,poi_1,poi_2,poi_3,poi_4,poi_5);
        Act_on_Prism_best(Description,Act,x_face_type,y_face_type,
                          poi_0,poi_1,poi_2,poi_3,poi_4,poi_5);
      }
    }
  }


  //  stop = false;
  //  stop = true;
  //  Boundary cell : 2.1.
  // ----------------------
  // Zweibein
  // Bemerkung: Test ob Kanten im Gebiet nicht notwendig wegen corners
  for(i=0;i<12 && !stop;++i) {
    for(k=0;k<2;++k) {
      // Punkt am Eck des Zweibeins
      corner = Transform((Edges_cell)i,(dir1D)k);

      poi_D.edge_point = corner_poi_typ;
      poi_D.corner = corner;

      // Richtungen des Zweibeins:
      dirA = opposite3D(OrthoA((Edges_cell)i));
      dirB = opposite3D(OrthoB((Edges_cell)i));
      dirC = Trans_dir((Edges_cell)i,(dir1D)k);

      // Zwei Ecken des Zweibeins
      poi_oppA.edge_point = poi_oppB.edge_point = corner_poi_typ;

      oppA = Next_corner(corner,dirA);
      oppB = Next_corner(corner,dirB);

      poi_oppA.corner = oppA;
      poi_oppB.corner = oppB;

      // Kantenpunkte des Zweibeins
      poi_A.edge_point  = poi_B.edge_point  = edge_poi_typ;
      poi_A2.edge_point = poi_B2.edge_point = edge_poi_typ;
      poi_C.edge_point  = edge_poi_typ;

      poi_A.corner  = oppA; poi_A.d  = dirB;
      poi_A2.corner = oppA; poi_A2.d = opposite3D(dirC);
      poi_B.corner  = oppB; poi_B.d  = dirA;
      poi_B2.corner = oppB; poi_B2.d = opposite3D(dirC);
      poi_C.corner  = corner; poi_C.d  = opposite3D(dirC);

      if(Description->edge_point(poi_A) && Description->edge_point(poi_A2) &&
         Description->edge_point(poi_B) && Description->edge_point(poi_B2) &&
         Description->edge_point(poi_C) &&
         corners[corner] && corners[oppA] && corners[oppB] ) {
        corners[corner] = corners[oppA] = corners[oppB] = false;
        //      stop = true;
        ori=0;
        // Orientierung:
        ori = 0;
        if(dirC==Wdir) ori=0;
        if(dirC==Edir) ori=1;
        if(dirC==Ndir) ori=0;
        if(dirC==Sdir) ori=1;
        if(dirC==Tdir) ori=1;
        if(dirC==Ddir) ori=0;
        
        if(((OneDpart(dirA)+OneDpart(dirB))&1)==ori) {
          poi_5 = poi_D;
          poi_2 = poi_C;
          
          poi_4 = poi_oppA;
          poi_6 = poi_oppB;

          poi_8 = poi_A;          poi_1 = poi_A2;
          poi_7 = poi_B;          poi_3 = poi_B2;       
        }
        else {
          poi_5 = poi_D;
          poi_2 = poi_C;

          poi_6 = poi_oppA;
          poi_4 = poi_oppB;

          poi_7 = poi_A;          poi_3 = poi_A2;
          poi_8 = poi_B;          poi_1 = poi_B2;
        }
  
        dirA = direction_from_to(poi_4.corner,poi_5.corner);
        dirB = poi_1.d;
        A_face_type = Calc_f_type(dirA,dirB);

        dirA = direction_from_to(poi_6.corner,poi_5.corner);
        dirB = poi_3.d;
        B_face_type = Calc_f_type(dirA,dirB);

        // Case 1.
        if(parallel_version)
          Tell_about_round_of_problem(Description->h(poi_1.corner, poi_1.d),
                                      Description->h(poi_8.corner, poi_8.d));
        if(Description->h(poi_1.corner, poi_1.d) <
           Description->h(poi_8.corner, poi_8.d) ||
           A_face_type == spitz) {
          // Anstelle: Act->Pyramid(poi_8,poi_5,poi_2,poi_1,poi_4); 
          if(A_face_type==an) {
            Act->Tetraeder(poi_8,poi_5,poi_1,poi_4);
            Act->Tetraeder(poi_8,poi_5,poi_2,poi_1);
          }
          else {
            Act->Tetraeder(poi_8,poi_4,poi_5,poi_2);
            Act->Tetraeder(poi_8,poi_4,poi_2,poi_1);
          }
          // Case 1.1.
          if(parallel_version)
            Tell_about_round_of_problem(Description->h(poi_7.corner, poi_7.d),
                                        Description->h(poi_3.corner, poi_3.d));
          if(Description->h(poi_7.corner, poi_7.d) <
             Description->h(poi_3.corner, poi_3.d) &&
             B_face_type == an) {
            Act->Tetraeder(poi_7,poi_8,poi_5,poi_3);
            Act->Tetraeder(poi_6,poi_7,poi_5,poi_3);
            Act->Tetraeder(poi_5,poi_3,poi_8,poi_2);
          }
          else {
            // Case 1.2.
            Act->Tetraeder(poi_5,poi_7,poi_8,poi_2);
            // Anstelle: Act->Pyramid(poi_7,poi_6,poi_3,poi_2,poi_5);
            if(B_face_type==an) {
              Act->Tetraeder(poi_7,poi_6,poi_3,poi_5);
              Act->Tetraeder(poi_7,poi_5,poi_3,poi_2);
            }
            else {
              Act->Tetraeder(poi_7,poi_6,poi_3,poi_2);
              Act->Tetraeder(poi_7,poi_6,poi_2,poi_5);
            }
          }
        }
        else {
          // Case 2.
          // Case 2.1.
          if(parallel_version)
            Tell_about_round_of_problem(Description->h(poi_3.corner, poi_3.d),
                                        Description->h(poi_7.corner, poi_7.d));
          if(Description->h(poi_3.corner, poi_3.d) <
             Description->h(poi_7.corner, poi_7.d) ||
             B_face_type == spitz) {
            Act->Tetraeder(poi_5,poi_1,poi_2,poi_7);
            Act->Tetraeder(poi_5,poi_7,poi_8,poi_1);
            Act->Tetraeder(poi_5,poi_4,poi_1,poi_8);

            // Anstelle: Act->Pyramid(poi_7,poi_5,poi_2,poi_3,poi_6);
            if(B_face_type==an) {
              Act->Tetraeder(poi_7,poi_6,poi_3,poi_5);
              Act->Tetraeder(poi_7,poi_5,poi_3,poi_2);
            }
            else {
              Act->Tetraeder(poi_7,poi_6,poi_3,poi_2);
              Act->Tetraeder(poi_7,poi_6,poi_2,poi_5);
            }
          }
          else {
            // Case 2.2.
            Act->Tetraeder(poi_5,poi_1,poi_2,poi_3);
            Act->Tetraeder(poi_5,poi_6,poi_7,poi_3);
            Act->Tetraeder(poi_5,poi_7,poi_8,poi_3);
            Act->Tetraeder(poi_5,poi_8,poi_1,poi_3);
            Act->Tetraeder(poi_5,poi_4,poi_1,poi_8);
          }
        }
      }
    }
  }
}



#endif
   

---