ug4
assembled_linear_operator_impl.h
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2011-2015: G-CSC, Goethe University Frankfurt
3  * Author: Andreas Vogel
4  *
5  * This file is part of UG4.
6  *
7  * UG4 is free software: you can redistribute it and/or modify it under the
8  * terms of the GNU Lesser General Public License version 3 (as published by the
9  * Free Software Foundation) with the following additional attribution
10  * requirements (according to LGPL/GPL v3 §7):
11  *
12  * (1) The following notice must be displayed in the Appropriate Legal Notices
13  * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14  *
15  * (2) The following notice must be displayed at a prominent place in the
16  * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17  *
18  * (3) The following bibliography is recommended for citation and must be
19  * preserved in all covered files:
20  * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21  * parallel geometric multigrid solver on hierarchically distributed grids.
22  * Computing and visualization in science 16, 4 (2013), 151-164"
23  * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24  * flexible software system for simulating pde based models on high performance
25  * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26  *
27  * This program is distributed in the hope that it will be useful,
28  * but WITHOUT ANY WARRANTY; without even the implied warranty of
29  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30  * GNU Lesser General Public License for more details.
31  */
32 
33 #ifndef __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__ASSEMBLED_LINEAR_OPERATOR_IMPL__
34 #define __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__ASSEMBLED_LINEAR_OPERATOR_IMPL__
35 
38 
39 #define PROFILE_ASS
40 #ifdef PROFILE_ASS
41  #define ASS_PROFILE_FUNC() PROFILE_FUNC()
42  #define ASS_PROFILE_BEGIN(name) PROFILE_BEGIN_GROUP(name, "discretization")
43  #define ASS_PROFILE_END() PROFILE_END()
44 #else
45  #define ASS_PROFILE_FUNC()
46  #define ASS_PROFILE_BEGIN(name)
47  #define ASS_PROFILE_END()
48 #endif
49 
50 namespace ug{
51 
52 template <typename TAlgebra>
53 void
55 {
56  if(m_spAss.invalid())
57  UG_THROW("AssembledLinearOperator: Assembling routine not set.");
58 
59 // assemble matrix (depending on u, i.e. J(u))
60  try{
61  m_spAss->assemble_jacobian(*this, u, m_gridLevel);
62  }
63  UG_CATCH_THROW("AssembledLinearOperator: Cannot assemble Jacobi matrix.");
64 }
65 
66 // Initialize the operator
67 template <typename TAlgebra>
68 void
70 {
71  if(m_spAss.invalid())
72  UG_THROW("AssembledLinearOperator: Assembling routine not set.");
73 
74 // create vector dummy
75  vector_type dummy;
76 
77 // assemble only matrix
78  try{
79  m_spAss->assemble_linear(*this, dummy, m_gridLevel);
80  }
81  UG_CATCH_THROW("AssembledLinearOperator::init: Cannot assemble Matrix.");
82 }
83 
84 // Initialize the operator
85 template <typename TAlgebra>
86 void
88 {
89 // todo: check that assembling is linear
90 
91  if(m_spAss.invalid())
92  UG_THROW("AssembledLinearOperator: Assembling routine not set.");
93 
94 // assemble matrix and rhs in one loop
95  try{
96  m_spAss->assemble_linear(*this, b, m_gridLevel);
97  }
98  UG_CATCH_THROW("AssembledLinearOperator::init_op_and_rhs:"
99  " Cannot assemble Matrix and Rhs.");
100 }
101 
102 template <typename TAlgebra>
103 void
105 {
106 #ifdef UG_PARALLEL
107  if(!c.has_storage_type(PST_CONSISTENT))
108  UG_THROW("Inadequate storage format of Vector c.");
109 #endif
110 
111 // perform check of sizes
112  if(c.size() != this->num_cols() || d.size() != this->num_rows())
113  UG_THROW("AssembledLinearOperator::apply: Size of matrix A ["<<
114  this->num_rows() << " x " << this->num_cols() << "] must match the "
115  "sizes of vectors x ["<<c.size()<<"], b ["<<d.size()<<"] for the "
116  " operation b = A*x. Maybe the operator is not initialized ?");
117 
118 // Apply Matrix
119  base_type::apply(d, c);
120 }
121 
122 // Compute d := d - J(u)*c
123 template <typename TAlgebra>
124 void
126 {
127 #ifdef UG_PARALLEL
128  if(!d.has_storage_type(PST_ADDITIVE))
129  UG_THROW("Inadequate storage format of Vector d.");
130  if(!c.has_storage_type(PST_CONSISTENT))
131  UG_THROW("Inadequate storage format of Vector c.");
132 #endif
133 
134 // check sizes
135  if(c.size() != this->num_cols() || d.size() != this->num_rows())
136  UG_THROW("AssembledLinearOperator::apply_sub: Size of matrix A ["<<
137  this->num_rows() << " x " << this->num_cols() << "] must match the "
138  "sizes of vectors x ["<<c.size()<<"], b ["<<d.size()<<"] for the "
139  " operation b -= A*x. Maybe the operator is not initialized ?");
140 
141 // Apply Matrix
142  base_type::matmul_minus(d,c);
143 }
144 
145 
146 template <typename TAlgebra>
148 {
149 // checks
150  if(m_spAss.invalid())
151  UG_THROW("AssembledLinearOperator: Assembling routine not set.");
152 
153 // set dirichlet values etc.
154  try{
155  m_spAss->adjust_solution(u, m_gridLevel);
156  }
157  UG_CATCH_THROW("AssembledLinearOperator::set_dirichlet_values:"
158  " Cannot assemble solution.");
159 }
160 
163 
165 template <typename TAlgebra>
168  typename TAlgebra::vector_type& u,
169  typename TAlgebra::vector_type& b)
170 {
171  ASS_PROFILE_BEGIN(ASS_AssembleLinearOperatorRhsAndSolution);
172 
173 // initialize operator
174  ASS_PROFILE_BEGIN(ASS_InitOperatorAndRhs);
175  try{
176  op.init_op_and_rhs(b);
177  }UG_CATCH_THROW("Cannot init the operator (assembling failed).");
178  ASS_PROFILE_END();
179 
180 // sets the dirichlet values in the solution
181  ASS_PROFILE_BEGIN(ASS_SetDirValues);
182  try{
183  op.set_dirichlet_values(u);
184  } UG_CATCH_THROW("Cannot set the dirichlet values in the solution.");
185  ASS_PROFILE_END();
186 
187  ASS_PROFILE_END();
188 }
189 
190 
191 } // end namespace ug
192 
193 #endif /* __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__ASSEMBLED_LINEAR_OPERATOR_IMPL__ */
#define ASS_PROFILE_BEGIN(name)
Definition: assembled_linear_operator_impl.h:42
#define ASS_PROFILE_END()
Definition: assembled_linear_operator_impl.h:43
matrix operator based on the assembling of a problem
Definition: assembled_linear_operator.h:60
void init_op_and_rhs(vector_type &b)
initializes the operator and assembles the passed rhs vector
Definition: assembled_linear_operator_impl.h:87
TAlgebra::vector_type vector_type
Type of Vector.
Definition: assembled_linear_operator.h:66
virtual void apply_sub(vector_type &d, const vector_type &c)
Compute d := d - J(u)*c.
Definition: assembled_linear_operator_impl.h:125
virtual void init()
initialize the operator
Definition: assembled_linear_operator_impl.h:69
virtual void apply(vector_type &d, const vector_type &c)
compute d = J(u)*c (here, J(u) is a Matrix)
Definition: assembled_linear_operator_impl.h:104
void set_dirichlet_values(vector_type &u)
Set Dirichlet values.
Definition: assembled_linear_operator_impl.h:147
@ PST_CONSISTENT
Definition: parallel_storage_type.h:68
@ PST_ADDITIVE
Definition: parallel_storage_type.h:69
#define UG_CATCH_THROW(msg)
Definition: error.h:64
#define UG_THROW(msg)
Definition: error.h:57
the ug namespace
void AssembleLinearOperatorRhsAndSolution(AssembledLinearOperator< TAlgebra > &op, typename TAlgebra::vector_type &u, typename TAlgebra::vector_type &b)
help function to assemble a linear operator
Definition: assembled_linear_operator_impl.h:167
size_t num_rows() const
Definition: sparsematrix_interface.h:38
size_t num_cols() const
Definition: sparsematrix_interface.h:39