app-kersten/Eddy2D_E.hh

Go to the documentation of this file.

#ifndef Eddy2D_hh
#define Eddy2D_hh

#include "space.hh"
#include "basics/debug.hh"

using concepts::Real;
using concepts::Real2d;
using concepts::Cmplx;

//************************************************************ Standardwerte **
// Polynomgrad
uint p = 0;
// Anzahl gleichmaessiger Verfeinerungen
uint l = 0;
// Anzahl geometrischer Verfeinerung
uint g = 0;
// leitendes Objekt
bool leitend = true;
// Grunddatei des Gitters
std::string gitter = "E2"; // E2, Es, 4Q oder 2R
// soll geloest werden
bool solve = true;
// sollen die Matrizen gespeichert werden
bool matrixoutput = false;
// soll eine graphische Ausgabe erfolgen
uint graphik = 2;
// Stringstreams fuer Ausgabedateinamen
std::ostringstream sl, sp, sg;
// Standardpfad
std::string path = SOURCEDIR "/app-kersten/Wirbelstrom2D/";
// Dateiname fuer Matrizenausgabe
std::string filename;

void input(int& argc, char** argv, std::string& probname) {
  int opt;
  filename = path + probname + "_Matrices.m";
  while ((opt = getopt(argc, argv, "-l:p:g:L:n:s:mG:hH")) != EOF)
    switch(opt) {
    case 'l': l = atoi(optarg); break;
    case 'p': p = atoi(optarg); break;
    case 'g': g = atoi(optarg); break;
    case 'L': leitend = std::atoi(optarg)>0; break;
    case 'n': gitter = optarg; break;
    case 's': solve = std::atoi(optarg)>0; break;
    case 'm': matrixoutput = true; break;
    case 'G': {
      graphik = std::atoi(optarg);
      if (graphik < 2) graphik = 2;
      break;
    }
    default:
      if ((opt != 'h') && (opt != 'H'))
  std::cout << "Option '" << opt << "' unkown." << std::endl;
      std::cout << "Usage: " << argv[0] 
    << " [-l LEVEL] [-p POLY] [-g GEOMLEVEL]"
    << "[-L CONDUCT] [-n MESH_NAME] [-s SOLVE]"
    << "[-m]" << std::endl
    << "where" << std::endl
    << "  LEVEL: level of refinement" << std::endl
    << "  POLY: starting polynomial degree" << std::endl
    << "  GEOMLEVEL: additonal level of geometrical refinement"
    << std::endl
    << "  CONDUCT: (0) non conducting material (1) conducting" 
    << std::endl
    << "  MESH_NAME: name of the mesh (thickLshape)" << std::endl
    << "  SOLVE: (0) not solving (1) solving system" << std::endl
    << "-m : matrices are stored into " << filename << std::endl;
  exit(1);
  break;
      }
  sl << l; sp << p; sg << g;
  std::cout << "Polynomgrad  : " << p << std::endl
      << "Verfeinerung : " << l << std::endl
            << "geom.Verfein : " << g << std::endl
            << "leitend      : " << leitend << std::endl
      << "Gitter       : " << gitter << std::endl
      << "Matrixausgabe: ";
  if (matrixoutput)
    std::cout << filename << std::endl;
  else {
    std::cout << "nein" << std::endl;
    filename.clear();
  }
}

//****************************************************************** Gitter **

std::string coord, elem, bound; // Gitterdateien
std::auto_ptr<concepts::Mesh2> msh;     // Gitter

void Gitter() {
  coord = path + "Coord_" + gitter + ".dat";
  elem  = path + "Elem_" + gitter + ".dat";
  bound = path + "Bound_" + gitter + ".dat";

  msh.reset
    (new concepts::Import2dMesh(coord.c_str(), elem.c_str(), bound.c_str()));

  std::cout << "Gitter aus den Dateien: " << std::endl
      << "  " << coord << std::endl
      << "  " << elem << std::endl
      << "  " << bound << std::endl;

  std::cout << "Gitter: " << *msh << std::endl;
}

//*************************************************************** Konstanten **

// Kreisfrequenz
const Real omega = 50.0/2/M_PI;   // Kreisfrequenz in 1/s

// Materialkonstanten
const Real sigma = 5.8e7;         // Leitfähigkeit in 1/Ohm/m
const Real eps_0 = 8.85419e-12;   // Dielektrizitätskonstante in s/(Ohm*m)
// const Real eps_0 = 1;
const Real mu_0  = 1.25644e-6;    // Permabilitätskonstante in Ohm*s/m
const Real J0 = 1;                // Stromdichte in A/m^2

//********************************************************** Materialformeln **

// stückweise definierte und konstante Leitfähigkeit
// Standardwert: 0
concepts::PiecewiseConstFormula<Real> Sigma;
// stückweise definierte Dielektrizität
// Standardwert: eps_0 (überall)
concepts::PiecewiseConstFormula<Real> Eps(eps_0);
// stückweise definierte und konstanter Quellstrom
// Standardwert: (0,0)
const Real sqrt2 = std::sqrt(0.5);
concepts::PiecewiseConstFormula<Real> J0x, J0y;

void material() {
  if (leitend) {
    Sigma[2] = sigma;                 // Werkstueck
    //     for(uint i = 2; i < 11; ++i) Sigma[i] = sigma;
  }
  std::cout << "Sigma: " << Sigma << std::endl;

  J0x[3] = J0x[4] = J0y[4] = J0y[5] = J0y[6] = J0x[10] = J0;
  J0y[3] = J0x[5] = J0x[9] = J0y[7] = 0.0;
  J0x[6] = J0x[7] = J0x[8] = J0y[8] = J0y[9] = J0y[10] = -J0; 
// J0x[3] = J0; J0y[3] = 0.0;
// J0x[4] = J0y[4] = sqrt2 * J0;
// J0x[5] = 0.0; J0y[5] = J0;
// J0x[6] = -sqrt2 * J0; J0y[6] = sqrt2 * J0;
// J0x[7] = -J0x[3]; J0x[8] = -J0x[4]; J0x[9] = -J0x[5]; J0x[10] = -J0x[6];
// J0y[7] = -J0y[3]; J0y[8] = -J0y[4]; J0y[9] = -J0y[5]; J0y[10] = -J0y[6];
  std::cout << "J0x: " << J0x << std::endl;
  std::cout << "J0y: " << J0y << std::endl;
}

//**************************************************************** FEM-Raum **

// Gemischter Raum
vectorial::Space<Real> Spc(2);

// geometrisches Verfeinern zu Ecken oder Kanten mit Attribut attrib
template<class F>
void verfeinern(F& spc, uint attrib, int pMax[2]) {
  for(uint i=0; i < g; ++i) {
    hp2D::APrioriRefinement refineSpace(spc, attrib, attrib, pMax);
    concepts::GlobalPostprocess<Real> post(spc);
    post(refineSpace);
    spc.rebuild();
    std::cout << std::endl;
  }
}

template<class F>
void verfeinerElm(F& spc, uint nrElm, uint l, uint p, bool rebuild = false) {
  concepts::AdaptiveAdjustP<2> a(l, p);
  std::auto_ptr<concepts::Space<Real>::Scanner> sc(spc.scan());
  const concepts::Element<Real>* elm;
  while (*sc && nrElm-- > 0)
    elm = &(*sc)++;
  spc.adjust(*elm, a);
  if (rebuild) spc.rebuild();
}

void elemente(concepts::Space<Real>& spc){
  std::auto_ptr<concepts::Space<Real>::Scanner> sc(spc.scan());
  uint i = 1;
  while (*sc) {
    const concepts::Element<Real>& elm = (*sc)++;
    std::cout << i++ << ".Element = " << elm << std::endl;
  }
}

//*********************************************************** Bilinearformen **

// a(.,.) Bilinearform
hp2Dedge::RotRot rotrot;
// m_Eddy(.,.) Bilinearform
std::auto_ptr<hp2Dedge::Identity> eddy;

// Systemmatrix
std::auto_ptr<concepts::SparseMatrix<Cmplx> > S;
// Rot-Rot-Matrix
std::auto_ptr<concepts::SparseMatrix<Real> > A;
// Wirbelstrom-Massenmatrix
std::auto_ptr<concepts::SparseMatrix<Real> > M_eddy;

std::auto_ptr<concepts::Vector<Cmplx> > rhs, sol, solN;

void rotrotMatrix(hp2Dedge::Space& spc) {
  concepts::ResourceMonitor timer;
  A.reset(new concepts::SparseMatrix<Real>(spc, rotrot));
  *A *= 1.0/mu_0;
  std::cout << "time for RotRot= " << timer << std::endl;
  std::cout << "A = " << *A << std::endl;
  A->storeMatlab(filename.c_str(), "A");
  A->addInto(*S, Cmplx(1.0,0.0));
}

void eddyMatrix(hp2Dedge::Space& spc) {
  // m_Eddy(.,.) Bilinearform mit Leitfaehigkeit
  eddy.reset(new hp2Dedge::Identity(Sigma));
  concepts::ResourceMonitor timer;
  M_eddy.reset(new concepts::SparseMatrix<Real>(spc, *eddy));
  *M_eddy *= omega;
  std::cout << "time for M_eddy= " << timer << std::endl;
  std::cout << "M_eddy = " << *M_eddy << std::endl;
  M_eddy->storeMatlab(filename.c_str(), "M_eddy", true);
  M_eddy->addInto(*S, Cmplx(0.0, 1.0));
}

void regMatrix(hp2Dedge::Space& spcE, hp2D::Space& spcN) {
  // Regularisierung
  vectorial::BilinearForm<Real> Bf(2);
  hp2Dedge::Graduv graduv;
  Bf.put(graduv, 1, 0);
  // Massenmatrix fuer Regularisierung
  concepts::ResourceMonitor timer;
  concepts::SparseMatrix<Real> M(Spc, Bf);
  std::cout << "time for Graduv= " << timer << std::endl;
  std::cout << "M = " << M << std::endl;
  M.storeMatlab(filename.c_str(), "M1", true);
  M.addInto(*S, Cmplx(10000000.0,0.0), 0, 0);
  M.addIntoT(*S, Cmplx(10000000.0,0.0), 0, 0); // transposed 
}

void lumpingMatrix(hp2Dedge::Space& spcE, hp2D::Space& spcN) {
  // Lumping-Beitrag
  hp2D::Laplace<Real> stiff;
  concepts::SparseMatrix<Real> A(spcN, stiff);
  std::cout << "A_ = " << A << std::endl;
  A.storeMatlab(filename.c_str(), "A_", true);
  A.addInto(*S, Cmplx(-1.0,0.0), spcE.dim(), spcE.dim());
}

void gesamtMatrix() {
  std::cout << "S = " << *S << std::endl;
  S->storeMatlab(filename.c_str(), "S", true);
//   S->compress();
//   std::cout << "S after compression = " << *S << std::endl;
}

template<class F>
void linearform(hp2Dedge::Space& spc, concepts::Space<F>& Spc) {
  rhs.reset(new concepts::Vector<Cmplx>(Spc));
  {
    hp2Dedge::Riesz lform(J0x, J0y);
    concepts::Vector<Real> rhs_(spc, lform);
    rhs_ *= -1.0;
    rhs->add(rhs_);
  }
  if (Spc.dim()<40)
    std::cout << "rhs = " << *rhs << std::endl;
}

void loeser(concepts::Space<Real>& spc) {
  sol.reset(new concepts::Vector<Cmplx>(spc));

  // kein Loesen, falls kein SuperLU vorhanden oder dies gewuenscht ist
#ifdef HAS_SuperLU
  if (!solve)
    return;

  concepts::SuperLU<Cmplx> solver(*S);
  solver(*rhs, *sol);

  if (spc.dim() < 40)
    std::cout << "sol = " << *sol << std::endl;
#endif
}

//******************************************************* graphische Ausgabe **

std::string dateiname(std::string varname, bool p = false) {
  std::string tmp(path);
  tmp += varname + '_' + gitter + "_l" + sl.str() + "_g" + sg.str();
  if (p) tmp += "_p" + sp.str();
  std::cout << tmp << std::endl;
  return tmp;
}

void formelausgabe(concepts::Space<Real>& spc,
                   concepts::PiecewiseFormulaBase<Real>& frm,
                   std::string frmstr) {
  graphics::MatlabGraphics(spc, dateiname(frmstr), frm);
}

void graphicalout(hp2Dedge::Space& spc, std::string probname) {
  if (graphik==0) return;

  hp2D::IntegrableQuad::rule().set(concepts::TRAPEZE, true, graphik);
  spc.recomputeShapefunctions();

  // Ausgabe von Gitter und des Materials in Matlab-Datei
  // Ansicht mit fill(x(msh),y(msh),u(msh))
  formelausgabe(spc, Sigma, "Sigma");
  formelausgabe(spc, J0x, "J0x");
  formelausgabe(spc, J0x, "J0y");

  graphics::drawMeshEPS(spc, 
      (dateiname(probname + "/Mesh/Mesh")+ ".eps").c_str());

  // Ausgabe der Lösung
  if (solve) {
    graphics::MatlabGraphics
      (spc, dateiname(probname + "/Grafik/" + probname + "_Sol", true), *sol);
  }
}

//**************************************************************** Energien **

// Verlustleistung 1/2 int j * e dx = 1/2 int sigma e * e dx [1 W]
template<class F>
Real dissipation(concepts::Vector<F>& sol) {
  concepts::Vector<F> tmp(sol.dim());
  (*M_eddy)(sol, tmp);
  tmp.apply(std::conj);
  return real(sol * tmp) / omega / 2.0;
}

// Magnetische Energie 1/2 int b * h dx [1 J]
template<class F>
Real magneticEnergy(concepts::Vector<F>& sol) {
  concepts::Vector<F> tmp(sol.dim());
  (*A)(sol, tmp);
  tmp.apply(std::conj);
  return real(sol * tmp) / 2.0 / (omega*omega);
}

#endif // Eddy2D_hh