v0.9.1
ConstrainMatrixCtx.hpp
Go to the documentation of this file.
1 /** \file ConstrainMatrixCtx.hpp
2  *
3  * Can be used if constrains are linear, i.e. are not function of time.
4  */
5 
6 /* This file is part of MoFEM.
7  * MoFEM is free software: you can redistribute it and/or modify it under
8  * the terms of the GNU Lesser General Public License as published by the
9  * Free Software Foundation, either version 3 of the License, or (at your
10  * option) any later version.
11  *
12  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
13  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  * License for more details.
16  *
17  * You should have received a copy of the GNU Lesser General Public
18  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
19 
20 #ifndef __PROJECTION_MATRIX_CTX_HPP__
21 #define __PROJECTION_MATRIX_CTX_HPP__
22 
23 /**
24  * \brief structure for projection matrices
25  * \ingroup projection_matrix
26  *
27  */
29 
31 
32  KSP kSP;
33  Mat C, CT, CCT, CTC, K;
35  Vec X, Qx, KQx;
37  bool createKSP;
41 
42  // Scatter is created form problem_x to problem_y, or scatter is given
43  // in the constructor
44 
45  VecScatter sCatter;
46  string xProblem, yProblem;
47 
48  PetscLogEvent MOFEM_EVENT_projInit;
49  PetscLogEvent MOFEM_EVENT_projQ;
50  PetscLogEvent MOFEM_EVENT_projP;
51  PetscLogEvent MOFEM_EVENT_projR;
52  PetscLogEvent MOFEM_EVENT_projRT;
53  PetscLogEvent MOFEM_EVENT_projCTC_QTKQ;
54 
55  /**
56  * Construct data structure to build operators for projection matrices
57  *
58  * User need to set matrix C to make it work
59  *
60  * \param x_problem problem on which vector is projected
61  * \param y_problem problem used to construct projection matrices
62  * \param create_ksp create ksp solver otherwise user need to set it up
63  */
64  ConstrainMatrixCtx(MoFEM::Interface &m_field, string x_problem,
65  string y_problem, bool create_ksp = true,
66  bool own_contrain_matrix = false);
67 
68  ConstrainMatrixCtx(MoFEM::Interface &m_field, VecScatter scatter,
69  bool create_ksp = true, bool own_contrain_matrix = false);
70 
71  virtual ~ConstrainMatrixCtx() {
72  ierr = destroyQorP();
73  CHKERRABORT(mField.get_comm(), ierr);
74  ierr = destroyQTKQ();
75  CHKERRABORT(mField.get_comm(), ierr);
76  if (ownConstrainMatrix) {
77  ierr = MatDestroy(&C);
78  CHKERRABORT(mField.get_comm(), ierr);
79  }
80  };
81 
82  PetscReal rTol, absTol, dTol;
83  PetscInt maxIts;
84 
85  /**
86  * \brief initialize vectors and matrices for Q and P shell matrices,
87  * scattering is set based on x_problem and y_problem
88  *
89  * \param x is a vector from problem x
90  */
92 
93  /**
94  * \brief initialize vectors and matrices for CTC+QTKQ shell matrices,
95  * scattering is set based on x_problem and y_problem
96  */
98 
99  /**
100  * \brief re-calculate CT and CCT if C matrix has been changed since
101  * initialization
102  */
104 
105  /**
106  * \brief re-calculate CTC matrix has been changed since initialization
107  */
109 
110  /**
111  * \brief destroy sub-matrices used for shell matrices P, Q, R, RT
112  */
114 
115  /**
116  * \brief destroy sub-matrices used for shell matrix QTKQ
117  */
119 
120  friend MoFEMErrorCode ProjectionMatrixMultOpQ(Mat Q, Vec x, Vec f);
121  friend MoFEMErrorCode ConstrainMatrixMultOpP(Mat P, Vec x, Vec f);
122  friend MoFEMErrorCode ConstrainMatrixMultOpR(Mat R, Vec x, Vec f);
123  friend MoFEMErrorCode ConstrainMatrixMultOpRT(Mat RT, Vec x, Vec f);
124  friend MoFEMErrorCode ConstrainMatrixMultOpCTC_QTKQ(Mat CTC_QTKQ, Vec x,
125  Vec f);
126 
129 };
130 
131 /**
132  * \brief Multiplication operator for Q = I-CTC(CCT)^-1C
133  *
134  * \code
135  * Mat Q; //for problem
136  * ConstrainMatrixCtx
137  projection_matrix_ctx(m_field,problem_name,contrains_problem_name);
138  * CHKERR MatCreateShell(PETSC_COMM_WORLD,m,m,M,M,&projection_matrix_ctx,&Q);
139  * CHKERR
140  MatShellSetOperation(Q,MATOP_MULT,(void(*)(void))ProjectionMatrixMultOpQ);
141  * CHKERR
142  MatShellSetOperation(Q,MATOP_DESTROY,(void(*)(void))ConstrainMatrixDestroyOpPorQ);
143  *
144  * \endcode
145 
146  * \ingroup projection_matrix
147 
148  */
149 MoFEMErrorCode ProjectionMatrixMultOpQ(Mat Q, Vec x, Vec f);
150 
151 /**
152  * \brief Multiplication operator for P = CT(CCT)^-1C
153  *
154  * \code
155  * Mat P; //for problem
156  * ConstrainMatrixCtx
157  projection_matrix_ctx(m_field,problem_name,contrains_problem_name);
158  * CHKERR MatCreateShell(PETSC_COMM_WORLD,m,m,M,M,&projection_matrix_ctx,&P);
159  * CHKERR
160  MatShellSetOperation(P,MATOP_MULT,(void(*)(void))ConstrainMatrixMultOpP);
161  *
162  * \endcode
163 
164  * \ingroup projection_matrix
165 
166  */
167 MoFEMErrorCode ConstrainMatrixMultOpP(Mat P, Vec x, Vec f);
168 
169 /**
170  * \brief Multiplication operator for R = CT(CCT)^-1
171  *
172  * \code
173  * Mat R; //for problem
174  * ConstrainMatrixCtx
175  projection_matrix_ctx(m_field,problem_name,contrains_problem_name);
176  * CHKERR MatCreateShell(PETSC_COMM_WORLD,m,m,M,M,&projection_matrix_ctx,&R);
177  * CHKERR
178  MatShellSetOperation(R,MATOP_MULT,(void(*)(void))ConstrainMatrixMultOpR);
179  *
180  * \endcode
181 
182  * \ingroup projection_matrix
183 
184  */
185 MoFEMErrorCode ConstrainMatrixMultOpR(Mat R, Vec x, Vec f);
186 
187 /**
188  * \brief Multiplication operator for RT = (CCT)^-TC
189  *
190  * \code
191  * Mat RT; //for problem
192  * ConstrainMatrixCtx
193  projection_matrix_ctx(m_field,problem_name,contrains_problem_name);
194  * CHKERR MatCreateShell(PETSC_COMM_WORLD,m,m,M,M,&projection_matrix_ctx,&RT);
195  * CHKERR
196  MatShellSetOperation(RT,MATOP_MULT,(void(*)(void))ConstrainMatrixMultOpRT);
197  *
198  * \endcode
199 
200  * \ingroup projection_matrix
201 
202  */
203 MoFEMErrorCode ConstrainMatrixMultOpRT(Mat RT, Vec x, Vec f);
204 
205 /**
206  * \brief Multiplication operator for RT = (CCT)^-TC
207  *
208  * \code
209  * Mat CTC_QTKQ; //for problem
210  * ConstrainMatrixCtx
211  projection_matrix_ctx(m_field,problem_name,contrains_problem_name);
212  * CHKERR
213  MatCreateShell(PETSC_COMM_WORLD,m,m,M,M,&projection_matrix_ctx,&CTC_QTKQ);
214  * CHKERR
215  MatShellSetOperation(CTC_QTKQ,MATOP_MULT,(void(*)(void))ConstrainMatrixMultOpCTC_QTKQ);
216  * CHKERR
217  MatShellSetOperation(CTC_QTKQ,MATOP_DESTROY,(void(*)(void))ConstrainMatrixDestroyOpQTKQ);
218  *
219  * \endcode
220  *
221 
222  * \ingroup projection_matrix
223 
224  */
225 MoFEMErrorCode ConstrainMatrixMultOpCTC_QTKQ(Mat CTC_QTKQ, Vec x, Vec f);
226 
227 /**
228  * \brief Destroy shell matrix Q
229  *
230  * \code
231  * Mat Q; //for problem
232  * ConstrainMatrixCtx
233  projection_matrix_ctx(m_field,problem_name,contrains_problem_name);
234  * CHKERR MatCreateShell(PETSC_COMM_WORLD,m,m,M,M,&projection_matrix_ctx,&Q);
235  * CHKERR
236  MatShellSetOperation(Q,MATOP_MULT,(void(*)(void))ProjectionMatrixMultOpQ);
237  * CHKERR
238  MatShellSetOperation(Q,MATOP_DESTROY,(void(*)(void))ConstrainMatrixDestroyOpPorQ);
239  *
240  * \endcode
241 
242  * \ingroup projection_matrix
243 
244  */
246 
247 /**
248  * \brief Destroy shell matrix
249  *
250  * \code
251  * Mat CTC_QTKQ; //for problem
252  * ConstrainMatrixCtx
253  projection_matrix_ctx(m_field,problem_name,contrains_problem_name);
254  * CHKERR MatCreateShell(PETSC_COMM_WORLD,m,m,M,M,&projection_matrix_ctx,&Q);
255  * CHKERR
256  MatShellSetOperation(Q,MATOP_MULT,(void(*)(void))ConstrainMatrixMultOpCTC_QTKQ);
257  * CHKERR
258  MatShellSetOperation(Q,MATOP_DESTROY,(void(*)(void))mat_destroy_QTKQ);
259  *
260  * \endcode
261 
262  * \ingroup projection_matrix
263 
264  */
266 
267 #endif // __PROJECTION_MATRIX_CTX_HPP__
268 
269 /**
270  \defgroup projection_matrix Constrain Projection Matrix
271  \ingroup user_modules
272 **/
PetscLogEvent MOFEM_EVENT_projRT
Deprecated interface functions.
friend MoFEMErrorCode ProjectionMatrixMultOpQ(Mat Q, Vec x, Vec f)
Multiplication operator for Q = I-CTC(CCT)^-1C.
friend MoFEMErrorCode ConstrainMatrixMultOpP(Mat P, Vec x, Vec f)
Multiplication operator for P = CT(CCT)^-1C.
PetscLogEvent MOFEM_EVENT_projQ
MoFEMErrorCode destroyQorP()
destroy sub-matrices used for shell matrices P, Q, R, RT
MoFEM::Interface & mField
structure for projection matrices
MoFEMErrorCode ConstrainMatrixDestroyOpQTKQ(Mat QTKQ)
Destroy shell matrix.
PetscLogEvent MOFEM_EVENT_projR
friend MoFEMErrorCode ConstrainMatrixMultOpCTC_QTKQ(Mat CTC_QTKQ, Vec x, Vec f)
Multiplication operator for RT = (CCT)^-TC.
MoFEMErrorCode recalculateCTandCCT()
re-calculate CT and CCT if C matrix has been changed since initialization
MoFEMErrorCode destroyQTKQ()
destroy sub-matrices used for shell matrix QTKQ
MoFEMErrorCode ConstrainMatrixMultOpCTC_QTKQ(Mat CTC_QTKQ, Vec x, Vec f)
Multiplication operator for RT = (CCT)^-TC.
MoFEMErrorCode ConstrainMatrixMultOpP(Mat P, Vec x, Vec f)
Multiplication operator for P = CT(CCT)^-1C.
friend MoFEMErrorCode ConstrainMatrixDestroyOpPorQ()
MoFEMErrorCode ConstrainMatrixDestroyOpPorQ(Mat Q)
Destroy shell matrix Q.
friend MoFEMErrorCode ConstrainMatrixMultOpRT(Mat RT, Vec x, Vec f)
Multiplication operator for RT = (CCT)^-TC.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
PetscLogEvent MOFEM_EVENT_projP
MoFEMErrorCode initializeQorP(Vec x)
initialize vectors and matrices for Q and P shell matrices, scattering is set based on x_problem and ...
MoFEMErrorCode initializeQTKQ()
initialize vectors and matrices for CTC+QTKQ shell matrices, scattering is set based on x_problem and...
MoFEMErrorCode ConstrainMatrixMultOpR(Mat R, Vec x, Vec f)
Multiplication operator for R = CT(CCT)^-1.
MoFEMErrorCode ProjectionMatrixMultOpQ(Mat Q, Vec x, Vec f)
Multiplication operator for Q = I-CTC(CCT)^-1C.
PetscLogEvent MOFEM_EVENT_projInit
MoFEMErrorCode ConstrainMatrixMultOpRT(Mat RT, Vec x, Vec f)
Multiplication operator for RT = (CCT)^-TC.
friend MoFEMErrorCode ConstrainMatrixMultOpR(Mat R, Vec x, Vec f)
Multiplication operator for R = CT(CCT)^-1.
virtual MPI_Comm & get_comm() const =0
MoFEMErrorCode recalculateCTC()
re-calculate CTC matrix has been changed since initialization
friend MoFEMErrorCode ConstrainMatrixDestroyOpQTKQ()
ConstrainMatrixCtx(MoFEM::Interface &m_field, string x_problem, string y_problem, bool create_ksp=true, bool own_contrain_matrix=false)
PetscLogEvent MOFEM_EVENT_projCTC_QTKQ