Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 : #include "libmesh/libmesh_config.h"
19 : #ifdef LIBMESH_HAVE_PETSC
20 :
21 : // Local includes
22 : #include "libmesh/petsc_shell_matrix.h"
23 : #include "libmesh/petsc_matrix_shell_matrix.h"
24 :
25 : namespace libMesh
26 : {
27 : template <typename T>
28 264 : PetscShellMatrix<T>::~PetscShellMatrix()
29 : {
30 132 : this->clear();
31 264 : }
32 :
33 : template <typename T>
34 0 : void PetscShellMatrix<T>::vector_mult (NumericVector<T> & dest,
35 : const NumericVector<T> & arg) const
36 : {
37 0 : PetscVector<T> & petsc_dest = cast_ref<PetscVector<T> &>(dest);
38 0 : const PetscVector<T> & petsc_arg = cast_ref<const PetscVector<T> &>(arg);
39 :
40 0 : LibmeshPetscCall(MatMult(_mat, petsc_arg.vec(), petsc_dest.vec()));
41 0 : }
42 :
43 :
44 :
45 : template <typename T>
46 0 : void PetscShellMatrix<T>::vector_mult_add (NumericVector<T> & dest,
47 : const NumericVector<T> & arg) const
48 : {
49 0 : PetscVector<T> & petsc_dest = cast_ref<PetscVector<T> &>(dest);
50 0 : const PetscVector<T> & petsc_arg = cast_ref<const PetscVector<T> &>(arg);
51 :
52 0 : LibmeshPetscCall(MatMultAdd(_mat, petsc_arg.vec(), petsc_dest.vec(), petsc_dest.vec()));
53 0 : }
54 :
55 :
56 : template <typename T>
57 132 : void PetscShellMatrix<T>::clear ()
58 : {
59 132 : if (this->initialized())
60 : {
61 : // If we encounter an error here, print a warning but otherwise
62 : // keep going since we may be recovering from an exception.
63 132 : PetscErrorCode ierr = MatDestroy (&_mat);
64 : if (ierr)
65 : libmesh_warning("Warning: MatDestroy returned a non-zero error code which we ignored.");
66 :
67 132 : this->_is_initialized = false;
68 : }
69 132 : }
70 :
71 : template <typename Obj>
72 622 : void init_shell_mat(Obj & obj,
73 : const numeric_index_type m,
74 : const numeric_index_type n,
75 : const numeric_index_type m_l,
76 : const numeric_index_type n_l,
77 : const numeric_index_type blocksize_in)
78 : {
79 : // Clear initialized matrices
80 622 : if (obj.initialized())
81 0 : obj.clear();
82 :
83 622 : PetscInt m_global = static_cast<PetscInt>(m);
84 622 : PetscInt n_global = static_cast<PetscInt>(n);
85 622 : PetscInt m_local = static_cast<PetscInt>(m_l);
86 622 : PetscInt n_local = static_cast<PetscInt>(n_l);
87 622 : PetscInt blocksize = static_cast<PetscInt>(blocksize_in);
88 :
89 622 : LibmeshPetscCall2(obj.comm(), MatCreate(obj.comm().get(), &obj._mat));
90 622 : LibmeshPetscCall2(obj.comm(), MatSetSizes(obj._mat, m_local, n_local, m_global, n_global));
91 622 : LibmeshPetscCall2(obj.comm(), MatSetBlockSize(obj._mat, blocksize));
92 622 : LibmeshPetscCall2(obj.comm(), MatSetType(obj._mat, MATSHELL));
93 :
94 : // Is prefix information available somewhere? Perhaps pass in the system name?
95 622 : LibmeshPetscCall2(obj.comm(), MatSetOptionsPrefix(obj._mat, ""));
96 622 : LibmeshPetscCall2(obj.comm(), MatSetFromOptions(obj._mat));
97 622 : LibmeshPetscCall2(obj.comm(), MatSetUp(obj._mat));
98 622 : LibmeshPetscCall2(obj.comm(), MatShellSetContext(obj._mat, &obj));
99 :
100 622 : obj._is_initialized = true;
101 622 : }
102 :
103 : template <typename Obj>
104 622 : void init_shell_mat(Obj & obj)
105 : {
106 14 : libmesh_assert(obj._dof_map);
107 :
108 622 : numeric_index_type my_m = obj._dof_map->n_dofs();
109 14 : numeric_index_type m_l = obj._dof_map->n_local_dofs();
110 622 : if (obj._omit_constrained_dofs)
111 : {
112 132 : my_m -= obj._dof_map->n_constrained_dofs();
113 132 : m_l -= obj._dof_map->n_local_constrained_dofs();
114 : }
115 :
116 14 : const numeric_index_type my_n = my_m;
117 14 : const numeric_index_type n_l = m_l;
118 :
119 622 : init_shell_mat(obj, my_m, my_n, m_l, n_l, obj._dof_map->block_size());
120 622 : }
121 :
122 : template <typename T>
123 132 : void PetscShellMatrix<T>::init ()
124 : {
125 132 : init_shell_mat(*this);
126 132 : }
127 :
128 : template <typename T>
129 264 : bool PetscShellMatrix<T>::initialized() const
130 : {
131 264 : return _is_initialized;
132 : }
133 :
134 : template <typename T>
135 264 : Mat PetscShellMatrix<T>::mat()
136 : {
137 264 : libmesh_error_msg_if(!_mat, "A petsc shell matrix is not created yet. Please call init() first.");
138 264 : return _mat;
139 : }
140 :
141 : //------------------------------------------------------------------
142 : // Explicit instantiations
143 : template class LIBMESH_EXPORT PetscShellMatrix<Number>;
144 : template void init_shell_mat(PetscShellMatrix<Number> & obj);
145 : template void init_shell_mat(PetscMatrixShellMatrix<Number> & obj);
146 : template void init_shell_mat(PetscShellMatrix<Number> & obj,
147 : const numeric_index_type m,
148 : const numeric_index_type n,
149 : const numeric_index_type m_l,
150 : const numeric_index_type n_l,
151 : const numeric_index_type blocksize_in);
152 : template void init_shell_mat(PetscMatrixShellMatrix<Number> & obj,
153 : const numeric_index_type m,
154 : const numeric_index_type n,
155 : const numeric_index_type m_l,
156 : const numeric_index_type n_l,
157 : const numeric_index_type blocksize_in);
158 :
159 : } // namespace libMesh
160 :
161 : #endif
|
