Assemble.hpp 15.9 KB
Newer Older
1
2
3
#ifndef _ASSEMBLE_H_
#define _ASSEMBLE_H_

4
5
#include "Algebra.hpp"
#include "Transfer.hpp"
6
#include "Parallel.hpp"
7
#include "Boundary.hpp"
8
#include "TimeSeries.hpp"
9

10
11
12
13
14
15
16
17
/*
 *  Interface defining the assembly of the global system of equations arising from local finite element equations.
 *  Most functions are declared twice: one global variant which is used by the Newton solver and one local variant
 *  calculating the contribution of the corresponding cell.
 *
 *  This interface should only be used for static problems. For time-dependent problems, the ITimeAssemble interface
 *  should be used.
 */
18

19
20
class IAssemble {
protected:
21
  int verbose = 1;
22
public:
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
  virtual const char *Name() const = 0;

  /*
   * Initializes the solution Vector u by setting the boundary conditions.
   */
  virtual void Initialize(Vector &u) const = 0;

  /*
   *  Computes the energy functional of a given vector u.
   */
  virtual double Energy(const Vector &u) const {
    double energy = 0;
    for (cell c = u.cells(); c != u.cells_end(); ++c) {
      Energy(c, u, energy);
    }
    return sqrt(PPM->Sum(energy));
  };

  virtual void Energy(const cell &c, const Vector &u, double &energy) const {
    Exit("void Energy(const cell &c, const Vector &u, double &energy)"
         " not implemented in your assembling.")
  };

  /*
   *  Assembles the first derivative of the energy
   *  function with respect to the test function v_j.
   *  The resulting residual is used by the Newton
   *  method as an error indicator.
   *  Todo: Check Comment and naming
   */
  virtual double Residual(const Vector &u, Vector &defect) const {
    defect = 0;
    for (cell ce = u.cells(); ce != u.cells_end(); ++ce)
      Residual(ce, u, defect);
    // TODO: Was mach ClearDirichletValues?
    defect.ClearDirichletValues();
    // TODO: Stattdessen PPM->Sum();?
    //MakeAdditive(defect);
61
    defect.Collect();
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
    return defect.norm();
  }

  virtual void Residual(const cell &c, const Vector &u, Vector &defect) const {
    Exit("void Residual(const cell &c, const Vector &u, Vector &defect)"
         " not implemented in your assembling.")
  };

  /*
   *  Assembles the second derivative of the energy
   *  function with respect to test functions v_j, w_j.
   *  The Jacobi matrix corresponds to the one used in the Newton method.
   */
  virtual void Jacobi(const Vector &u, Matrix &jacobi) const {
    jacobi = 0;
    for (cell c = jacobi.cells(); c != jacobi.cells_end(); ++c)
      Jacobi(c, u, jacobi);
    jacobi.ClearDirichletValues();
  }

  virtual void Jacobi(const cell &c, const Vector &u, Matrix &jacobi) const {
    Exit("void Jacobi(const cell &c, const Vector &u, Matrix &jacobi)"
         " not implemented in your assembling.")
  };

  /*
   * Can be used to print information about your assembled system.
   */
  virtual void PrintInfo() const {
91
92
93
    mout.PrintInfo("Assemble", verbose,
                   PrintInfoEntry("Name", this->Name(), 0)
                   );
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
  }

  virtual void PrintInfo(const Vector &u) const {
    Exit("void PrintInfo(const Vector &u) not implemented in your assembling.")
  }

  /*
   * Defines a projection of the fine mesh onto the coarse mesh
   */
  virtual void Project(const Vector &fineVector, Vector &coarseVector) const {}

  /*
   * TODO: Somebody who uses this should document it
   */
  virtual void AssembleTransfer(TransferMatrix &) const {}
109
110
111
112
113
114
};

/*
 * Todo move to own file
 */
class MatrixTransfer : public Transfer {
115
116
117
  const IAssemble &assemble;
  TransferMatrixGraph *TG;
  TransferMatrix *TM;
118
public:
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
  MatrixTransfer(const IAssemble &A)
      : assemble(A), TG(0), TM(0) {}

  void multiply(Vector &f, const Vector &c) const {
    vout(5) << "AssembleTransfer: multiply" << endl;
    vout(6) << "Prolongate coarse " << endl << c;
    f = *TM * c;
    f.ClearDirichletValues();
    vout(7) << "Prolongate fine " << endl << f;
  }

  void multiply_transpose(Vector &c, const Vector &f) const {
    vout(5) << "AssembleTransfer: multiply_transpose" << endl;
    vout(7) << "Restrict fine " << endl << f;
    c = f * (*TM);
    c.ClearDirichletValues();
135
    c.Collect();
136
137
138
139
140
141
142
    vout(6) << "Restrict coarse " << endl << c;
  }

  void Project(const Vector &f, Vector &c) const {
    vout(5) << "AssembleTransfer: Project" << endl;
    vout(7) << "Project fine " << endl << f;
    assemble.Project(f, c);
143
    c.Average();
144
145
146
    vout(6) << "Project coarse " << endl << c;
  }

147
  void Construct(const VectorMatrixBase &fg, const VectorMatrixBase &cg) {
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
    if (TG == 0) {
      TG = new TransferMatrixGraph(fg, cg);
      TM = new TransferMatrix(*TG);
      assemble.AssembleTransfer(*TM);
    }
  }

  void Destruct() {
    if (TG) {
      delete TG;
      TG = 0;
      delete TM;
      TM = 0;
    }
  }

  Transfer *GetTransferPointer() const {
    return new MatrixTransfer(assemble);
  }

  friend Logging &operator<<(Logging &s, const MatrixTransfer &T) {
    return s << *(T.TG);
  }
171
172
};

173
class ILinearTimeAssemble {
174
175
protected:
  int verbose = 1;
176

177
  int step = 0;
178

179
  double t_ = 0.0;
180

181
  double dt_ = 0.0;
182

183
184
  TimeSeries timeSeries;

185
public:
186
  virtual const char *Name() const = 0;
187

188
  virtual void Initialize(Vector &u) = 0;
189

190
191
192
193
194
  virtual void RHS(double t, Vector &rhs) const = 0;

  virtual void MassMatrix(Matrix &massMatrix) const = 0;

  virtual void SystemMatrix(Matrix &systemMatrix) const = 0;
195

196
  virtual void PrintInfo(const Vector &u) const {
197
    // Todo print info about vector at current time step
198
  };
199

200
201
202
203
    virtual void PrintInfo() const {
        // Todo print info about time series for example
    };

204
  TimeSeries &GetTimeSeries() { return timeSeries; }
205

206
  ILinearTimeAssemble(TimeSeries timeSeries) : timeSeries(timeSeries) {}
207

208
  virtual ~ILinearTimeAssemble() = default;
209

210
211
212
213
  virtual void InitTimeStep() {
    dt_ = -t_;
    t_ = timeSeries.NextTStep(t_);
    dt_ += t_;
214
  }
215

216
  virtual void FinishTimeStep(double t, Vector &u_new) {
217
218
    ++step;
  }
219

220
221
  int Step() {
    return step;
222
  }
223

224
225
226
227
228
229
230
  double Time() {
    return t_;
  }

  double TimeStep() {
    return dt_;
  }
231

232
  virtual double Energy(const Vector &u) const {
233
    Exit("Energy not implemented")
234
235
  }
};
236

237
238
239
240
class TAssemble : public IAssemble {
  bool symmetric;
public:
  const char *Name() const override { return "TAssemble"; }
241

242
  const Vector &U_old() const { return *u_; }
243

244
  const Vector &V_old() const { return *v_; }
245

246
  const Vector &A_old() const { return *a_; }
247

248
249
250
251
252
253
254
255
256
257
258
protected:
  const Vector *u_ = nullptr;
  const Vector *v_ = nullptr;
  const Vector *a_ = nullptr;

  int step = 0;
  double beta = 0.0;
  double gamma = 0.0;
  double t_ = 0.0;
  double dt_ = 0.0;
public:
259

260
261
262
263
264
265
266
267
268
269
270
  void Initialize(Vector &u) const override {
    Dirichlet(t_, u);
  }

  void Dirichlet(Vector &u) const {
    Dirichlet(t_, u);
  }

  virtual void Dirichlet(double t, Vector &u) const {
    u.ClearDirichletFlags();
    for (cell c = u.cells(); c != u.cells_end(); ++c) Dirichlet(t, c, u);
271
    u.DirichletConsistent();
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
  }

  virtual void Dirichlet(double t, const cell &c, Vector &u) const {
    Exit("Dirichlet not implemented")
  }

  double Residual(const Vector &u, Vector &r) const override {
    return Residual(t_, u, r);
  }

  virtual double Residual(double t, const Vector &u, Vector &r) const {
    r = 0;
    for (cell c = u.cells(); c != u.cells_end(); ++c)
      Residual(t, c, u, r);
    r.ClearDirichletValues();
287
    r.Collect();
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
    return r.norm();
  }

  virtual void Residual(double t, const cell &c, const Vector &u, Vector &b) const {
    Exit("Residual not implemented")
  }

  void Jacobi(const Vector &u, Matrix &jacobi) const override {
    Jacobi(t_, u, jacobi);
  }

  virtual void Jacobi(double t, const Vector &u, Matrix &A) const {
    A = 0;
    for (cell c = u.cells(); c != u.cells_end(); ++c) Jacobi(t, c, u, A);
    if (symmetric) A.Symmetric();
    A.ClearDirichletValues();
  }

  virtual void Jacobi(double t, const cell &c, const Vector &u, Matrix &A) const {
    Exit("Jacobi not implemented")
  }

  virtual void PrintInfo(const Vector &u) const {};

  explicit TAssemble(bool sym = false) : symmetric(sym) {}

  void resetTAssemble() {
    u_ = nullptr;
    v_ = nullptr;
    a_ = nullptr;

    step = 0;
    beta = 0.0;
    gamma = 0.0;
    t_ = 0.0;
    dt_ = 0.0;
  }

  int getStep() {
    return step;
  }

  virtual ~TAssemble() = default;

  virtual void Initialize(double t) {
    t_ = t;
  }

  virtual void Initialize(double t, Vector &u) {
    t_ = t;
  }

  virtual void Initialize(double t, const Vector &u, const Vector &v, Vector &a) {
    t_ = t;
  }

  virtual void BeforeInitTimeStep(double t_old, double t, double t_new,
                                  const Vector &u_old, const Vector &u,
                                  Vector &u_new) {}

  virtual void InitTimeStep(double t, double dt) {
    t_ = t;
    dt_ = dt;
  }

  virtual void InitTimeStep(double t_old, double t, double t_new,
                            const Vector &u_old, const Vector &u,
                            Vector &u_new) {
    BeforeInitTimeStep(t_old, t, t_new, u_old, u, u_new);
    t_ = t_new;
    dt_ = t_new - t;
    u_ = &u;
  }

  virtual void InitTimeStep(double t_old, double t, double t_new,
                            const Vector &u_old, Vector &u,
                            Vector &u_new) {
    BeforeInitTimeStep(t_old, t, t_new, u_old, u, u_new);
    t_ = t_new;
    dt_ = t_new - t;
    u_ = &u;
  }

  virtual void InitTimeStep(double t_old, double t, double t_new,
                            const Vector &u_old, const Vector &u,
                            const Vector &u_new, const Vector &v_old,
                            const Vector &a_old) {
    t_ = t_new;
    dt_ = t_new - t;
    u_ = &u_old;
    v_ = &v_old;
    a_ = &a_old;
  }

  virtual void BeforeFinishTimeStep(double t, const Vector &u_new,
                                    bool SpecialTimeStep = false) {}

  virtual void AfterFinishTimeStep(double t, const Vector &u_new,
                                   bool SpecialTimeStep = false) {}

  virtual void FinishTimeStep(double t, const Vector &u_new,
                              bool SpecialTimeStep = false) {
    BeforeFinishTimeStep(t, u_new, SpecialTimeStep);
    ++step;
    AfterFinishTimeStep(t, u_new, SpecialTimeStep);
  }

  virtual void FinishTimeStep(double t, Vector &u_new, bool SpecialTimeStep = false) {
    ++step;
  }

  virtual void FinishTimeStep(double t, const Vector &u_new, const Vector &v_new,
                              bool SpecialTimeStep = false) {
    ++step;
  }
403

404
405
406
  virtual double NextTimeStep(const Vector &u, double t, double t_new) const {
    return t_new;
  }
407

408
  virtual void Finalize(double t, const Vector &u_new) {}
409

410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
  virtual void Finalize(double t, const Vector &u_new,
                        const Vector &v_new,
                        const Vector &a_new) {}

  virtual bool BlowUp(const Vector &u_new) const { return false; }

  virtual void MassMatrix(Matrix &massMatrix) const {};

  virtual void MassMatrix(const Vector &u, Matrix &massMatrix) const {};

  virtual void SystemMatrix(Matrix &systemMatrix) const {};

  virtual void SystemMatrix(double t,
                            const Vector &u,
                            Matrix &systemMatrix) const {};
425

426
  virtual void RHS(Vector &rhs) const {};
427

428
  virtual void RHS(double t, Vector &rhs) const {};
429

430
  virtual void RHS(double t, const Vector &u, Vector &rhs) const {};
431

432
433
434
  virtual double Energy(double t, const cell &c, const Vector &u) const {
    Exit("double Energy(double t, const cell &c, const Vector &u) not implemented")
  };
435

436
437
438
439
440
441
442
  virtual double Energy(double t, const Vector &u) const {
    double E = 0;
    for (cell c = u.cells(); c != u.cells_end(); ++c) {
      E += Energy(t, c, u);
    }
    return PPM->Sum(E);
  }
443

444
445
446
  double Energy(const Vector &u) const override {
    return Energy(t_, u);
  }
447

448
449
450
451
452
453
  int Step() const { return step; }

  void SetNewmarkParameters(double _gamma, double _beta) {
    beta = _beta;
    gamma = _gamma;
  }
454
455
};

456
457
class INonLinearTimeAssemble : public TAssemble {};

458
class HAssemble : public IAssemble {
459
460
461
462
  double t_;
  int step;
  mutable int level;
  bool symmetric;
463
private:
464
  virtual const char *Name() const { return "HAssemble"; }
465

466
  virtual int find_level(const Vector &u) const { return 0; }
467

468
469
470
  void read_level(const Vector &u) const {
    level = find_level(u);
  }
471
472

public:
473
  HAssemble(bool sym = false) : IAssemble(), symmetric(sym), step(0), t_(0) {}
474

475
  virtual ~HAssemble() {}
476

477
  int Level() const { return level; }
478

479
480
481
482
  void Dirichlet(Vector &u) const {
    read_level(u);
    Dirichlet(t_, u);
  }
483

484
485
486
487
  double Energy(Vector &u) const {
    read_level(u);
    return Energy(t_, u);
  }
488

489
490
491
492
  double Residual(const Vector &u, Vector &b) const {
    read_level(u);
    return Residual(t_, u, b);
  }
493

494
495
496
497
498
499
  void Jacobi(const Vector &u, Matrix &A) const {
    mout.StartBlock("HAssemble");
    read_level(u);
    Jacobi(t_, u, A);
    mout.EndBlock(TimeLevel, 1);
  }
500

501
502
503
  void AfterNewtonDefect(const Vector &u) const {
    CountPlasticPoints(u);
  }
504

505
  virtual void NewtonMinusCorrection(Vector &u, const Vector &c) {}
506

507
  virtual void NewtonPlusCorrection(Vector &u, const Vector &c) {}
508

509
  virtual void NewtonUpdate(Vector &u, const Vector &c) {}
510

511
  virtual void NewtonBackdate(Vector &u, const Vector &c) {}
512

513
  virtual void ExtrapolateRotationVector(double s) {}
514

515
  virtual void ExtrapolateRotationVector(double s, const Vector &u) {}
516

517
  virtual void ResetRotationVector() {}
518

519
  virtual void ResetRotationVector(const Vector &u) {}
520

521
  virtual void project_eulerian_angle(Vector &c) {}
522
523

private:
524
525
526
527
528
529
  virtual void Dirichlet(double t, const cell &c, Vector &u) const {}

  virtual void Dirichlet(double t, Vector &u) const {
    u.ClearDirichletFlags();
    for (cell c = u.cells(); c != u.cells_end(); ++c)
      Dirichlet(t, c, u);
530
    u.DirichletConsistent();
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
  }

  virtual bool detailed_energies() const { return false; }

  virtual void reset_energies() const {}

  virtual double sum_energies() const { return 0; }

  virtual double Energy(double t, const cell &c, Vector &u) const { return 0.0; }

  virtual double Energy(double t, Vector &u) const {
    if (detailed_energies()) {
      reset_energies();
      for (cell c = u.cells(); c != u.cells_end(); ++c)
        Energy(t, c, u);
      return sum_energies();
    } else {
      double E = 0;
      for (cell c = u.cells(); c != u.cells_end(); ++c)
        E += Energy(t, c, u);
      E = PPM->Sum(E);
      return E;
    }
  }

  virtual void Residual(double t, const cell &c, const Vector &u,
                        Vector &b) const = 0;

  virtual double Residual(double t, const Vector &u, Vector &b) const {
    b = 0;
    for (cell c = u.cells(); c != u.cells_end(); ++c)
      Residual(t, c, u, b);
    b.ClearDirichletValues();
564
    b.Collect();
565
566
567
568
569
570
571
572
573
574
575
576
    return b.norm();
  }

  virtual void Jacobi(double t, const cell &c, const Vector &u,
                      Matrix &A) const = 0;

  virtual void Jacobi(double t, const Vector &u, Matrix &A) const {
    A = 0;
    for (cell c = u.cells(); c != u.cells_end(); ++c) Jacobi(t, c, u, A);
    if (symmetric) A.Symmetric();
    A.ClearDirichletValues();
  }
577
578

public:
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
  void Initialize(double t, Vector &u) {
    Init(t, u);
    t_ = t;
  }

  void InitTimeStep(double t_old, double t, double t_new,
                    const Vector &u_old, const Vector &u,
                    Vector &u_new) {
    InitTStep(t_old, t, t_new, u_old, u, u_new);
    t_ = t_new;
  }

  void Finalize(double t_new, const Vector &u_new) {
    Finish(t_new, u_new);
  }

  void FinishTimeStep(double t_new, const Vector &u_new,
                      bool SpecialTimeStep = true) {
    FinishTStep(t_new, u_new, SpecialTimeStep);
    ++step;
  }

  virtual void CountPlasticPoints(const Vector &u) const {}
602
603

private:
604
  virtual void Init(double t, Vector &u) {}
605

606
607
608
  virtual void InitTStep(double t_old, double t, double t_new,
                         const Vector &u_old, const Vector &u,
                         Vector &u_new) {}
609

610
  virtual void Finish(double t_new, const Vector &u_new) {}
611

612
613
  virtual void
  FinishTStep(double t_new, const Vector &u_new, bool SpecialTimeStep = true) {}
614
615

public:
616
  virtual bool BlowUp(const Vector &u_new) const { return false; }
617

618
  int Step() const { return step; }
619
620
621
};

#endif // of #ifndef _ASSEMBLE_H_