couple.cc 4.74 KB
Newer Older
1
2
3
4
5
6
7
8
#include "AMDiS.h"

using namespace std;
using namespace AMDiS;

class G : public AbstractFunction<double, WorldVector<double> >
{
public:
9
10
  double operator()(const WorldVector<double>& x) const 
  {
11
12
    return exp(-10.0 * (x * x));
  }
13
14
15
16
17
};

class F : public AbstractFunction<double, WorldVector<double> >
{
public:
18
  F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {}
19

20
21
  double operator()(const WorldVector<double>& x) const 
  {
22
    int dow = x.getSize();
23
24
25
26
    double r2 = (x * x);
    double ux = exp(-10.0 * r2);
    return -(400.0 * r2 - 20.0 * dow) * ux;
  }
27
28
29
30
31
32
33
34
};

class MyCoupledIteration : public ProblemIterationInterface
{
public:
  MyCoupledIteration(ProblemStatBase *prob1,
		     ProblemStatBase *prob2)
    : problem1(prob1),
35
36
37
      problem2(prob2),
      name("MyCoupledIteration")
  {}
38
39
40
41
42
43
44

  void beginIteration(AdaptInfo *adaptInfo) 
  {
    FUNCNAME("StandardProblemIteration::beginIteration()");
    MSG("\n");
    MSG("begin of iteration number: %d\n", adaptInfo->getSpaceIteration()+1);
    MSG("=============================\n");
45
  }
46
47
48
49
50
51

  void endIteration(AdaptInfo *adaptInfo) {
    FUNCNAME("StandardProblemIteration::endIteration()");
    MSG("\n");
    MSG("end of iteration number: %d\n", adaptInfo->getSpaceIteration()+1);
    MSG("=============================\n");
52
  }
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67

  Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo = FULL_ITERATION) 
  {
    Flag flag, markFlag;
    if(toDo.isSet(MARK)) markFlag = problem1->markElements(adaptInfo);
    if(toDo.isSet(ADAPT) && markFlag.isSet(MESH_REFINED)) flag = problem1->refineMesh(adaptInfo);

    if(toDo.isSet(BUILD)) problem1->buildAfterCoarsen(adaptInfo, markFlag);
    if(toDo.isSet(SOLVE)) problem1->solve(adaptInfo);
 
    if(toDo.isSet(BUILD)) problem2->buildAfterCoarsen(adaptInfo, markFlag);
    if(toDo.isSet(SOLVE)) problem2->solve(adaptInfo);

    if(toDo.isSet(ESTIMATE)) problem1->estimate(adaptInfo);    
    return flag;
68
  }
69
70
71
72

  int getNumProblems() 
  {
    return 2;
73
  }
74
75
76
77

  ProblemStatBase *getProblem(int number = 0) 
  {
    FUNCNAME("CoupledIteration::getProblem()");
78
79
80
81
82
83

    if (number == 0) 
      return problem1;
    if (number == 1) 
      return problem2;

84
85
    ERROR_EXIT("invalid problem number\n");
    return NULL;
86
87
88
89
90
91
92
93
94
95
96
97
  }

  const std::string& getName()
  {
    return name;
  }

  void serialize(std::ostream&)
  {}

  void deserialize(std::istream&)
  {}
98
99
100

private:
  ProblemStatBase *problem1;
101

102
  ProblemStatBase *problem2;
103
104

  std::string name;
105
106
107
108
109
};

class Identity : public AbstractFunction<double, double>
{
public:
110
  Identity(int degree) : AbstractFunction<double, double>(degree) {}
111

112
113
  double operator()(const double& x) const 
  {
114
115
    return x;
  }
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
};

int main(int argc, char* argv[])
{
  FUNCNAME("main");
  TEST_EXIT(argc == 2)("usage: couple initfile\n");
  Parameters::init(true, argv[1]);

  // ===== create and init the first problem ===== 
  ProblemScal problem1("problem1");
  problem1.initialize(INIT_ALL);

  // ===== create and init the second problem ===== 
  Flag initFlag = 
    INIT_FE_SPACE |
    INIT_SYSTEM |
    INIT_SOLVER |
    INIT_FILEWRITER;

  Flag adoptFlag =
    CREATE_MESH |
    INIT_MESH;

  ProblemScal problem2("problem2");
  problem2.initialize(initFlag,
		      &problem1,
		      adoptFlag);

  // ===== create operators for problem1 =====
145
  Operator matrixOperator1(Operator::MATRIX_OPERATOR, problem1.getFeSpace());
146
  matrixOperator1.addSecondOrderTerm(new Laplace_SOT);
147
148
  problem1.addMatrixOperator(&matrixOperator1);

149
150
  int degree = problem1.getFeSpace()->getBasisFcts()->getDegree();
  Operator rhsOperator1(Operator::VECTOR_OPERATOR, problem1.getFeSpace());
151
  rhsOperator1.addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree)));
152
153
154
  problem1.addVectorOperator(&rhsOperator1);

  // ===== create operators for problem2 =====
155
  Operator matrixOperator2(Operator::MATRIX_OPERATOR, problem2.getFeSpace());
156
  matrixOperator2.addZeroOrderTerm(new Simple_ZOT);
157
158
  problem2.addMatrixOperator(&matrixOperator2);
  
159
  Operator rhsOperator2(Operator::VECTOR_OPERATOR, problem2.getFeSpace());
160
161
  rhsOperator2.addZeroOrderTerm(new VecAtQP_ZOT(problem1.getSolution(), 
						new Identity(degree)));
162
163
  problem2.addVectorOperator(&rhsOperator2);

164
165
166
167
168
169
  // ===== add boundary conditions for problem1 =====
  problem1.addDirichletBC(1, new G);

  // ===== add boundary conditions for problem1 =====
  //problem2.addDirichletBC(1, new G);

170
  // ===== create adaptation loop and iteration interface =====
171
  AdaptInfo *adaptInfo = new AdaptInfo("couple->adapt", 1);
172
  MyCoupledIteration coupledIteration(&problem1, &problem2);
173
  AdaptStationary *adapt = new AdaptStationary("couple->adapt",
174
175
176
177
178
179
180
181
182
183
					       &coupledIteration,
					       adaptInfo);

  // ===== start adaptation loop =====
  adapt->adapt();

  // ===== write solution ===== 
  problem1.writeFiles(adaptInfo, true);
  problem2.writeFiles(adaptInfo, true);
}