src/expde/formulas/diffop.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>
#else
 #include <iostream>
 #include <fstream>
#endif 

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

// Include level 4:
#include "boundy.h"
#include "diffop.h"


Boundary_matrices::Boundary_matrices(int max_number_matrices,
                                     BoCeData* bocedata) {
  int i,j;
  Tetraeder_storage* tet;
  bool* needed;

  max_num_matrices = max_number_matrices;
  if(bocedata->Exists_bocellpoint())
    leng_col = bocedata->Give_number_points() + 1;
  else
    leng_col = bocedata->Give_number_points();
  leng_row = new int[leng_col];
  needed   = new bool[leng_col];

  // calc size of matrix and entires
  int max_leng;
  max_leng=0;
  for(i=0;i<leng_col;++i) {
    // a) size of row
    for(j=0;j<leng_col;++j) needed[j] = false;

    for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
      if(tet->N0()==i || tet->N1()==i || tet->N2()==i || tet->N3()==i) {
        needed[tet->N0()]=true;
        needed[tet->N1()]=true;
        needed[tet->N2()]=true;
        needed[tet->N3()]=true;
      }
    }
    leng_row[i] = 0;
    for(j=0;j<leng_col;++j)
      if(needed[j])
        ++leng_row[i];

    if(max_leng<leng_row[i]) max_leng=leng_row[i];
  }

  // b) entries
  num_row = new int[max_leng*leng_col];

  int num;
  for(i=0;i<leng_col;++i) {
    for(j=0;j<leng_col;++j) needed[j] = false;
    
    for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
      if(tet->N0()==i || tet->N1()==i || tet->N2()==i || tet->N3()==i) {
        needed[tet->N0()]=true;
        needed[tet->N1()]=true;
        needed[tet->N2()]=true;
        needed[tet->N3()]=true;
      }
    }
    num = 0;
    for(j=0;j<leng_col;++j)
      if(needed[j]) {
        ++num;
        num_row[i+leng_col*num] = j;
      }
  }
  
  // c) initialize matrix;
  matrix = new double[max_num_matrices*max_leng*leng_col];
}

---