EllipticPDESolver.cpp 2.41 KB
Newer Older
niklas.baumgarten's avatar
niklas.baumgarten committed
1
#include "EllipticPDESolver.hpp"
niklas.baumgarten's avatar
niklas.baumgarten committed
2
3


4
void EllipticPDESolver::run(SampleSolution &solution) {
niklas.baumgarten's avatar
niklas.baumgarten committed
5
  newton->operator()(*assemble, solution.U);
niklas.baumgarten's avatar
niklas.baumgarten committed
6
7
}

niklas.baumgarten's avatar
niklas.baumgarten committed
8
// Todo maybe multiply weights to Q
9
void EllipticPDESolver::computeQ(SampleSolution &solution) {
niklas.baumgarten's avatar
niklas.baumgarten committed
10
11
12
13
  if (quantity == "L2") solution.Q = assemble->L2(solution.U);
  else if (quantity == "H1") solution.Q = assemble->H1(solution.U);
  else if (quantity == "Energy") solution.Q = assemble->Energy(solution.U);
  else if (quantity == "Inflow") solution.Q = assemble->InflowOutflow(solution.U).first;
niklas.baumgarten's avatar
niklas.baumgarten committed
14
  else if (quantity == "Outflow") solution.Q = assemble->InflowOutflow(solution.U).second;
niklas.baumgarten's avatar
niklas.baumgarten committed
15
16
17
18
19
20
21
22
  else if (quantity == "L2Error") solution.Q = assemble->L2Error(solution.U);
  else if (quantity == "EnergyError") solution.Q = assemble->EnergyError(solution.U);
  else if (quantity == "L2CellAverageError")
    solution.Q = assemble->L2CellAverageError(solution.U);
  else if (quantity == "MaxError") solution.Q = assemble->MaxError(solution.U);
  else if (quantity == "FluxError") solution.Q = assemble->FluxError(solution.U);
  else if (quantity == "FaceError") solution.Q = assemble->FaceError(solution.U);
  else Exit("Quantity of interest not implemented")
niklas.baumgarten's avatar
niklas.baumgarten committed
23
24
25
}

void EllipticPDESolver::computeCost(SampleSolution &solution) {
niklas.baumgarten's avatar
niklas.baumgarten committed
26
  if (costMeasure == "size") solution.C = solution.U.size();
niklas.baumgarten's avatar
niklas.baumgarten committed
27
//    else if (costMeasure == "time") solution.Cost = solution.U.size(); // Todo
niklas.baumgarten's avatar
niklas.baumgarten committed
28
  else Exit("Cost measure not implemented")
niklas.baumgarten's avatar
niklas.baumgarten committed
29
30
31
}

void EllipticPDESolver::plotSolution(SampleSolution &solution) {
niklas.baumgarten's avatar
niklas.baumgarten committed
32
  if (!plotting) return;
33

niklas.baumgarten's avatar
niklas.baumgarten committed
34
35
  IDiscretization *pressureDisc;
  if (typeid(*assemble).name() == typeid(HybridEllipticAssemble).name())
niklas.baumgarten's avatar
niklas.baumgarten committed
36
    pressureDisc = new LagrangeDiscretization(GetDisc()->GetMeshes(), 0, 1);
niklas.baumgarten's avatar
niklas.baumgarten committed
37
  else
niklas.baumgarten's avatar
niklas.baumgarten committed
38
    pressureDisc = new LagrangeDiscretization(GetDisc()->GetMeshes(), 1, 1);
niklas.baumgarten's avatar
niklas.baumgarten committed
39

niklas.baumgarten's avatar
niklas.baumgarten committed
40
41
42
  SampleSolution pressure(pressureDisc, solution.Id(), "Pressure");
  assemble->SetPressure(solution.U, pressure.U);
  mpp::plot(pressure.IdString()) << pressure.U << mpp::endp;
43

niklas.baumgarten's avatar
niklas.baumgarten committed
44
  LagrangeDiscretization fluxDisc(GetDisc()->GetMeshes(), 0, SpaceDimension);
niklas.baumgarten's avatar
niklas.baumgarten committed
45
46
47
  SampleSolution flux(&fluxDisc, solution.Id(), "Flux");
  assemble->SetFlux(solution.U, flux.U);
  mpp::plot(flux.IdString()) << flux.U << mpp::endp;
niklas.baumgarten's avatar
niklas.baumgarten committed
48

niklas.baumgarten's avatar
niklas.baumgarten committed
49
  LagrangeDiscretization kappaDisc(GetDisc()->GetMeshes(), 0, 1);
50
51
52
  SampleSolution kappa(&kappaDisc, solution.Id(), "Kappa");
  assemble->GetProblem()->Permeability(kappa.U);
  mpp::plot(kappa.IdString()) << kappa.U << mpp::endp;
53

niklas.baumgarten's avatar
niklas.baumgarten committed
54
  delete pressureDisc;
niklas.baumgarten's avatar
niklas.baumgarten committed
55
}