coupled problem (polarization and navierStokes) added to extensions demos

extensions/demo/other/src/polarizationField_navierStokes.cc

+/******************************************************************************
+ *
+ * Extension of AMDiS - Adaptive multidimensional simulations
+ *
+ * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
+ * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
+ *
+ * Authors: Simon Praetorius et al.
+ *
+ * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+ * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+ *
+ *
+ * See also license.opensource.txt in the distribution.
+ * 
+ ******************************************************************************/
+
+#include "AMDiS.h"
+#include "PolarizationField.h"
+#include "NavierStokes_TaylorHood.h"
+
+#include "GenericOperatorTerm.h"
+
+// coupling structures
+#if __cplusplus > 199711L
+  #include "CouplingBaseProblem2_cxx11.h"
+#else
+  #include "CouplingBaseProblem2.h"
+#endif
+   
+using namespace AMDiS;
+
+
+///
+class MyPolarizationField : public PolarizationField
+{
+public:
+  typedef PolarizationField super;
+  
+public:
+  MyPolarizationField(std::string name)
+  : super(name) {}
+  
+  void solveInitialProblem(AdaptInfo *adaptInfo)
+  {
+    struct RandomNormalizedVector : public AbstractFunction<WorldVector<double>, WorldVector<double> >
+    {
+      RandomNormalizedVector() : AbstractFunction<WorldVector<double>, WorldVector<double> >(1) 
+      {
+	std::srand(time(0));
+      }
+      
+      WorldVector<double> operator()(const WorldVector<double>& x) const
+      {
+	WorldVector<double> p;
+	for (int i = 0; i < p.getSize(); i++)
+	  p[i] = cos(4.0*m_pi * ((std::rand() / static_cast<double>(RAND_MAX)) - 0.5));
+	
+	double nrm = norm(p);
+	for (int i = 0; i < p.getSize(); i++)
+	  p[i] *= 1.0/nrm;
+	return p;
+      }
+    };
+    
+    getVectorField()->interpol(new RandomNormalizedVector);
+    for (int i = 0; i < Global::getGeo(WORLD); i++)
+      transformDOF(getVectorField(), getSolution()->getDOFVector(i), new AMDiS::Component2<>(i));
+  }
+  
+  void fillBoundaryConditions() override
+  {
+    std::vector<DOFVector<double>*> values(3);
+    values[0] = new DOFVector<double>(getFeSpace(0), "v0"); values[0]->set(0.0);
+    values[1] = new DOFVector<double>(getFeSpace(0), "v1"); values[1]->set(1.0);
+    values[2] = new DOFVector<double>(getFeSpace(0), "v2"); values[2]->set(-1.0);
+    int idx = 1;
+    for (int i = 0; i < Global::getGeo(WORLD); i++) {
+      for (int j = 0; j < Global::getGeo(WORLD); j++) {
+	double value1 = double(i == j);
+	double value2 = -double(i == j);
+	prob->addDirichletBC(idx, 1-j + dow, 1-j, new AMDiS::Const<double, WorldVector<double> >(value1));
+	prob->addDirichletBC(idx+1, 1-j + dow, 1-j, new AMDiS::Const<double, WorldVector<double> >(value1));
+      }
+      idx += 2;
+    }
+  }
+};
+
+
+///
+struct MyNavierStokes : public ::detail::NavierStokes_TaylorHood<ExtendedProblemStat>
+{
+public:
+  typedef ::detail::NavierStokes_TaylorHood<ExtendedProblemStat> super;
+
+public:
+  MyNavierStokes(std::string name)
+    : super(name) {}
+    
+  void fillBoundaryConditions() override
+  {
+    struct DrivenCavityBC : AbstractFunction<double, WorldVector<double> >
+    {
+      double operator()(const WorldVector<double> &x) const 
+      {
+	return std::max(0.0, 1.0 - 4.0 * sqr(x[0] - 0.5));
+      }
+    };
+    
+    DOFVector<double> *zeroDOF = new DOFVector<double>(getFeSpace(0), "zero");
+    zeroDOF->set(0.0);
+
+    /// at rigid wall: no-slip boundary condition
+    for (size_t i = 0; i < dow; i++) {
+      prob->addDirichletBC(1, i, i, zeroDOF);
+      prob->addDirichletBC(2, i, i, zeroDOF);
+      prob->addDirichletBC(3, i, i, zeroDOF);
+    }
+
+    /// at upper wall: prescribed velocity
+    prob->addDirichletBC(4, 0, 0, new DrivenCavityBC);
+    prob->addDirichletBC(4, 1, 1, zeroDOF);
+  }
+  
+};
+
+
+///
+struct PolarNavierStokes : public CouplingBaseProblem<ProblemStat, MyPolarizationField, MyNavierStokes>
+{
+public: // typedefs
+  typedef CouplingBaseProblem<ProblemStat, MyPolarizationField, MyNavierStokes> super;
+  
+public: // methods
+  PolarNavierStokes(std::string name, MyPolarizationField& pProb_, MyNavierStokes& nsProb_)
+    : super(name, pProb_, nsProb_), 
+      pProb(pProb_), nsProb(nsProb_),
+      M0(1.0),
+      R0(1.0)
+  {
+    Parameters::get(name + "->M0", M0);
+    Parameters::get(name + "->R0", R0);
+  }
+    
+  void fillCouplingOperators() override
+  {
+    
+    WorldVector<DOFVector<double>* > P, PS;
+    for (size_t i = 0; i < dow; i++) {
+      P[i]= pProb.getSolution()->getDOFVector(i);
+      PS[i]= pProb.getSolution()->getDOFVector(i + dow);
+    }
+    
+    for (size_t i = 0; i < dow; i++) {
+      // < [(grad(P#)^T * P]_i , psi >
+      Operator *opPSGradP_i = new Operator(nsProb.getFeSpace(i), pProb.getFeSpace(0));
+      for (size_t j = 0; j < dow; j++)
+	addZOT(opPSGradP_i, M0 * valueOf(PS[j]) * derivativeOf(P[j], i) );
+      nsProb.getProblem(0)->addVectorOperator(opPSGradP_i, i);
+      
+      // < [div(P# o P - P o P#)]_i, psi > = < div(sigma^S)_i, psi >
+      Operator *opDivPPS = new Operator(nsProb.getFeSpace(i), pProb.getFeSpace(0));
+      addZOT(opDivPPS, ( 0.5*M0) * valueOf(P[1-i])  * derivativeOf(PS[i],   1-i) );
+      addZOT(opDivPPS, (-0.5*M0) * valueOf(P[i])    * derivativeOf(PS[1-i], 1-i) );
+      addZOT(opDivPPS, ( 0.5*M0) * valueOf(PS[i])   * derivativeOf(P[1-i],  1-i) );
+      addZOT(opDivPPS, (-0.5*M0) * valueOf(PS[1-i]) * derivativeOf(P[i],    1-i) );
+      nsProb.getProblem(0)->addVectorOperator(opDivPPS, i);
+    }
+    
+    
+    WorldVector<DOFVector<double>* > vel;
+    for (size_t i = 0; i < dow; i++)
+      vel[i]= nsProb.getSolution()->getDOFVector(i);
+    
+    for (size_t i = 0; i < dow; i++) {
+      // < u * grad(P_i) , psi >
+      for (size_t j = 0; j < dow; j++) {
+	Operator *opUGradP = new Operator(pProb.getFeSpace(0), pProb.getFeSpace(0));
+	addFOT(opUGradP, valueOf(vel[j]), j, GRD_PHI);
+	pProb.getProblem()->addMatrixOperator(opUGradP, i, i);
+      }
+      
+      // < -1/2 * [(grad(u) - grad(u)^T) * P]_i, psi >
+      // = < -S(grad(u),P)_i, psi >
+      Operator *opGradUP = new Operator(pProb.getFeSpace(0), pProb.getFeSpace(0));
+      addZOT(opGradUP, (R0 * (i == 0 ? 0.5 :-0.5)) * derivativeOf(vel[1], 0) );
+      addZOT(opGradUP, (R0 * (i == 0 ?-0.5 : 0.5)) * derivativeOf(vel[0], 1) );
+      pProb.getProblem()->addMatrixOperator(opGradUP, i, 1-i);
+    }
+    
+  }
+    
+protected:
+  MyPolarizationField& pProb; 
+  MyNavierStokes& nsProb;
+  
+  double M0;
+  double R0;
+};
+
+
+int main(int argc, char** argv)
+{ FUNCNAME("main");
+
+  AMDiS::init(argc, argv);  
+  Timer t;
+  
+  MyPolarizationField pProb("p");
+  MyNavierStokes nsProb("ns");
+  
+  PolarNavierStokes mainProb("main", pProb, nsProb);  
+  mainProb.initialize(INIT_ALL);
+
+  // Adapt-Infos
+  AdaptInfo adaptInfo("adapt", mainProb.getNumComponents());
+  AdaptInstationary adaptInstat("adapt", mainProb, adaptInfo, mainProb, adaptInfo);
+
+  mainProb.initTimeInterface();
+  int error_code = adaptInstat.adapt(); 
+
+  MSG("elapsed time= %d sec\n", t.elapsed());
+  AMDiS::finalize();
+
+  return error_code;
+}