src/expde/extemp/opoth.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
//

// ------------------------------------------------------------
// opoth.cc
//
// ------------------------------------------------------------

// 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 2:
#include "../basic/basic.h"

// Include level 3:
#include "../domain/domain.h"

// Include level 4:
#include "../formulas/boundy.h"
#include "../formulas/loc_sten.h"
#include "../formulas/diffop.h"
#include "../formulas/diffopc.h"
#include "../formulas/diffopv.h"

// Include level 5:
#include "../grid/gpar.h"
#include "../grid/parallel.h"
#include "../grid/mgcoeff.h"
#include "../grid/sto_man.h"
#include "../grid/gridbase.h"
#include "../grid/grid.h"
#include "../grid/input.h"


// Include level 5.1
#include "../evpar/evpar.h"


// Include level 6:
#include "variable.h"
#include "opera.h"
#include "operav.h"
#include "sim.h"
#include "extra.h"
#include "res_op.h"
#include "array.h"




               // 1. functions of sim.h
               // 2. functions of extra.h
               // 3. functions of array.h

               
// 1. functions of sim.h


DWrapSim<DExprSimLit>
operator== (Variable& a, double value)
{
  typedef DWrapSim<DExprSimLit> ExprT;
  return ExprT(DExprSimLit(a,value)); 
}



DWrapSim<DExprSimlocLit>
operator== (const Local_var& a, double x)
{
  typedef DWrapSim<DExprSimlocLit > ExprT;
  return ExprT(DExprSimlocLit(a,x)); 
}


DWrapSim<DExprSimloc<DExprVAR> >
operator== (const Local_var& a, Variable& b)
{
  typedef DWrapSim<DExprSimloc<DExprVAR> > ExprT;
  return ExprT(DExprSimloc<DExprVAR>(a,DExprVAR(b))); 
}


// 2. functions of extra.h

// Extrapolation at the boundary
DExpr<DVar_Extra_Op>
Extrapolate_near_boundary(Variable a)
{
  return DExpr<DVar_Extra_Op>(DVar_Extra_Op(a));
}



double DVar_Extra_Op::Give_at_Index(const Index3D I, const Grid* grid) const {
  int l;
  Index3D In, Inn, Innn;
  int i;
  Pointtype Ptypn, Ptypnn, Ptypnnn;
  double un, unn, unnn;

  l = grid->Max_level();

  //  return grid->Give_variable(I,l)[v_.Number_variable()];

  for(i=0;i<8;++i) {
    In   = I.next((dir_sons)i, l);
    Inn  = In.next((dir_sons)i, l);
    Innn = Inn.next((dir_sons)i, l);
    if(grid->Point_in_domain(In) && grid->Point_in_domain(Inn)
       && grid->Point_in_domain(Innn)) {
      Ptypn  = grid->Give_type(In);
      Ptypnn = grid->Give_type(Inn);
      Ptypnnn = grid->Give_type(Innn);
      if((interior == Ptypn   || uniform <= Ptypn) &&
         (interior == Ptypnn  || uniform <= Ptypnn) &&
         (interior == Ptypnnn || uniform <= Ptypnnn)) {
        un   = grid->Give_variable(In,l)[v_.Number_variable()];
        unn  = grid->Give_variable(Inn,l)[v_.Number_variable()];
        unnn = grid->Give_variable(Innn,l)[v_.Number_variable()];
        
        return unnn - 3.0*unn + 3.0 * un;
      }
    }
  }

  for(i=0;i<8;++i) {
    In  = I.next((dir_sons)i, l);
    Inn = In.next((dir_sons)i, l);
    if(grid->Point_in_domain(In) && grid->Point_in_domain(Inn)) {
      Ptypn  = grid->Give_type(In);
      Ptypnn = grid->Give_type(Inn);
      if((interior == Ptypn  || uniform <= Ptypn) &&
         (interior == Ptypnn || uniform <= Ptypnn)) {
        un  = grid->Give_variable(In,l)[v_.Number_variable()];
        unn = grid->Give_variable(Inn,l)[v_.Number_variable()];
        
        //      h_mesh = grid->H_mesh() / Zweipotenz(grid->Max_level());
        
        return 2.0*un - unn;
      }
    }
  }
  return grid->Give_variable(I,l)[v_.Number_variable()];

  /*

  weight = 0;
  sum    = 0.0;
  for(i=0;i<3;++i) for(j=0;j<3;++j) for(k=0;k<3;++k) {
    In = I.next((Ort1D)i, (Ort1D)j, (Ort1D)k, l);
    Ptyp = grid->Give_type(In);
    if(interior == Ptyp || uniform <= Ptyp) {
      weight = weight + 1;
      sum = sum + grid->Give_variable(In,l)[v_.Number_variable()];
    }
  }
  if(weight==0) {
    //    cout << "Null";
    return grid->Give_variable(I,l)[v_.Number_variable()];
  }
  //  cout << "Eins";
  return sum / weight;
  */
}

               
// 3. functions of array.h

Schleife For(Local_int& I, int begin, int end)  {
  return Schleife(I.Give_pointer(),begin,end);
}

SchleifeSum Sum(Local_int& I, int begin, int end)  {
  return SchleifeSum(I.Give_pointer(),begin,end);
}

---