MFEM  v3.3.2
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
operator.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 #include "vector.hpp"
13 #include "operator.hpp"
14 
15 #include <iostream>
16 #include <iomanip>
17 
18 namespace mfem
19 {
20 
21 void Operator::FormLinearSystem(const Array<int> &ess_tdof_list,
22  Vector &x, Vector &b,
23  Operator* &Aout, Vector &X, Vector &B,
24  int copy_interior)
25 {
26  const Operator *P = this->GetProlongation();
27  const Operator *R = this->GetRestriction();
28  Operator *rap;
29 
30  if (P)
31  {
32  // Variational restriction with P
33  B.SetSize(P->Width());
34  P->MultTranspose(b, B);
35  X.SetSize(R->Height());
36  R->Mult(x, X);
37  rap = new RAPOperator(*P, *this, *P);
38  }
39  else
40  {
41  // rap, X and B point to the same data as this, x and b
42  X.NewDataAndSize(x.GetData(), x.Size());
43  B.NewDataAndSize(b.GetData(), b.Size());
44  rap = this;
45  }
46 
47  if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
48 
49  // Impose the boundary conditions through a ConstrainedOperator, which owns
50  // the rap operator when P and R are non-trivial
51  ConstrainedOperator *A = new ConstrainedOperator(rap, ess_tdof_list,
52  rap != this);
53  A->EliminateRHS(X, B);
54  Aout = A;
55 }
56 
57 void Operator::RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
58 {
59  const Operator *P = this->GetProlongation();
60  if (P)
61  {
62  // Apply conforming prolongation
63  x.SetSize(P->Height());
64  P->Mult(X, x);
65  }
66  else
67  {
68  // X and x point to the same data
69  }
70 }
71 
72 void Operator::PrintMatlab(std::ostream & out, int n, int m) const
73 {
74  using namespace std;
75  if (n == 0) { n = width; }
76  if (m == 0) { m = height; }
77 
78  Vector x(n), y(m);
79  x = 0.0;
80 
81  out << setiosflags(ios::scientific | ios::showpos);
82  for (int i = 0; i < n; i++)
83  {
84  x(i) = 1.0;
85  Mult(x, y);
86  for (int j = 0; j < m; j++)
87  {
88  if (y(j))
89  {
90  out << j+1 << " " << i+1 << " " << y(j) << '\n';
91  }
92  }
93  x(i) = 0.0;
94  }
95 }
96 
97 
99  bool _own_A)
100  : Operator(A->Height(), A->Width()), A(A), own_A(_own_A)
101 {
102  constraint_list.MakeRef(list);
103  z.SetSize(height);
104  w.SetSize(height);
105 }
106 
108 {
109  w = 0.0;
110 
111  for (int i = 0; i < constraint_list.Size(); i++)
112  {
113  w(constraint_list[i]) = x(constraint_list[i]);
114  }
115 
116  A->Mult(w, z);
117 
118  b -= z;
119 
120  for (int i = 0; i < constraint_list.Size(); i++)
121  {
122  b(constraint_list[i]) = x(constraint_list[i]);
123  }
124 }
125 
126 void ConstrainedOperator::Mult(const Vector &x, Vector &y) const
127 {
128  if (constraint_list.Size() == 0)
129  {
130  A->Mult(x, y);
131  return;
132  }
133 
134  z = x;
135 
136  for (int i = 0; i < constraint_list.Size(); i++)
137  {
138  z(constraint_list[i]) = 0.0;
139  }
140 
141  A->Mult(z, y);
142 
143  for (int i = 0; i < constraint_list.Size(); i++)
144  {
145  y(constraint_list[i]) = x(constraint_list[i]);
146  }
147 }
148 
149 }
int Size() const
Logical size of the array.
Definition: array.hpp:110
void NewDataAndSize(double *d, int s)
Set the Vector data and size, deleting the old data, if owned.
Definition: vector.hpp:101
Array< int > constraint_list
List of constrained indices/dofs.
Definition: operator.hpp:387
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:320
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:42
virtual void Mult(const Vector &x, Vector &y) const
Constrained operator action.
Definition: operator.cpp:126
int Size() const
Returns the size of the vector.
Definition: vector.hpp:113
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
double * GetData() const
Definition: vector.hpp:121
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Definition: operator.hpp:52
virtual const Operator * GetProlongation() const
Prolongation operator from linear algebra (linear system) vectors, to input vectors for the operator...
Definition: operator.hpp:65
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:36
The operator x -&gt; R*A*P*x.
Definition: operator.hpp:320
virtual const Operator * GetRestriction() const
Restriction operator from input vectors for the operator to linear algebra (linear system) vectors...
Definition: operator.hpp:68
void SetSubVectorComplement(const Array< int > &dofs, const double val)
Set all vector entries NOT in the &#39;dofs&#39; array to the given &#39;val&#39;.
Definition: vector.cpp:601
Vector w
Auxiliary vectors.
Definition: operator.hpp:390
void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, Operator *&A, Vector &X, Vector &B, int copy_interior=0)
Form a constrained linear system using a matrix-free approach.
Definition: operator.cpp:21
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
Reconstruct a solution vector x (e.g. a GridFunction) from the solution X of a constrained linear sys...
Definition: operator.cpp:57
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:24
void EliminateRHS(const Vector &x, Vector &b) const
Eliminate &quot;essential boundary condition&quot; values specified in x from the given right-hand side b...
Definition: operator.cpp:107
void PrintMatlab(std::ostream &out, int n=0, int m=0) const
Prints operator with input size n and output size m in Matlab format.
Definition: operator.cpp:72
ConstrainedOperator(Operator *A, const Array< int > &list, bool own_A=false)
Constructor from a general Operator and a list of essential indices/dofs.
Definition: operator.cpp:98
void MakeRef(T *, int)
Make this Array a reference to a pointer.
Definition: array.hpp:644
Vector data type.
Definition: vector.hpp:41
Operator * A
The unconstrained Operator.
Definition: operator.hpp:388
Square Operator for imposing essential boundary conditions using only the action, Mult()...
Definition: operator.hpp:384
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:64
Abstract operator.
Definition: operator.hpp:21
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:25