Mercurial > fem-fenics-eugenio

comparison src/assemble.cc @ 247:8ca45824938e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add factories hierarchy for matrices' and vectors' assembly
author Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
date Sat, 12 Jul 2014 12:30:13 +0200
parents 66071811eef8
children 5e9b5bbdc56b
comparison
equal deleted inserted replaced
246:a2991e2a085b 247:8ca45824938e
15 this program; if not, see <http://www.gnu.org/licenses/>. 15 this program; if not, see <http://www.gnu.org/licenses/>.
16 */ 16 */
17 17
18 #include "form.h" 18 #include "form.h"
19 #include "boundarycondition.h" 19 #include "boundarycondition.h"
20 #include "femfenics_factory.h"
20 21
21 DEFUN_DLD (assemble, args, nargout, 22 DEFUN_DLD (assemble, args, nargout,
22 "-*- texinfo -*-\n\ 23 "-*- texinfo -*-\n\
23 @deftypefn {Function File} {[@var{A}], [@var{x}(Optional)]} = \ 24 @deftypefn {Function File} {[@var{A}], [@var{x}(Optional)]} = \
24 assemble (@var{form_a}, @var{DirichletBC}) \n\ 25 assemble (@var{form_a}, @var{DirichletBC}) \n\
77 const dolfin::Form & a = frm.get_form (); 78 const dolfin::Form & a = frm.get_form ();
78 a.check (); 79 a.check ();
79 80
80 if (a.rank () == 2) 81 if (a.rank () == 2)
81 { 82 {
82 dolfin::parameters["linear_algebra_backend"] = "uBLAS"; 83 femfenics_factory factory;
84
83 dolfin::Matrix A; 85 dolfin::Matrix A;
84 dolfin::assemble (A, a); 86 dolfin::assemble (A, a);
85 87
86 for (std::size_t i = 1; i < nargin; ++i) 88 for (std::size_t i = 1; i < nargin; ++i)
87 { 89 {
102 } 104 }
103 else 105 else
104 error ("assemble: unknown argument type"); 106 error ("assemble: unknown argument type");
105 } 107 }
106 108
107 // Get capacity of the dolfin sparse matrix 109 retval(0) = factory.matrix (A);
108 boost::tuples::tuple<const std::size_t*,
109 const std::size_t*,
110 const double*, int>
111 aa = A.data ();
112
113 int nnz = aa.get<3> ();
114 std::size_t nr = A.size (0), nc = A.size (1);
115 std::vector<double> data_tmp;
116 std::vector<std::size_t> cidx_tmp;
117
118 dim_vector dims (nnz, 1);
119 octave_idx_type nz = 0, ii = 0;
120 Array<octave_idx_type>
121 ridx (dims, 0),
122 cidx (dims, 0);
123 Array<double> data (dims, 0);
124
125 octave_idx_type* orow = ridx.fortran_vec ();
126 octave_idx_type* oc = cidx.fortran_vec ();
127 double* ov = data.fortran_vec ();
128
129 for (std::size_t i = 0; i < nr; ++i)
130 {
131 A.getrow (i, cidx_tmp, data_tmp);
132 nz += cidx_tmp.size ();
133
134 for (octave_idx_type j = 0;
135 j < cidx_tmp.size (); ++j)
136 {
137 orow [ii + j] = i;
138 oc [ii + j] = cidx_tmp [j];
139 ov [ii + j] = data_tmp [j];
140 }
141
142 ii = nz;
143 }
144
145 dims(0) = ii;
146 ridx.resize (dims);
147 cidx.resize (dims);
148 data.resize (dims);
149
150 SparseMatrix sm (data, ridx, cidx, nr, nc);
151 retval(0) = sm;
152 } 110 }
153 111
154 else if (a.rank () == 1) 112 else if (a.rank () == 1)
155 { 113 {
114 femfenics_factory factory;
115
156 dolfin::Vector A; 116 dolfin::Vector A;
157 dolfin::assemble (A, a); 117 dolfin::assemble (A, a);
158 118
159 for (std::size_t i = 1; i < nargin; ++i) 119 for (std::size_t i = 1; i < nargin; ++i)
160 { 120 {
174 } 134 }
175 else 135 else
176 error ("assemble: unknown argument type"); 136 error ("assemble: unknown argument type");
177 } 137 }
178 138
179 dim_vector dims; 139 retval(0) = factory.vector (A);
180 dims.resize (2);
181 dims(0) = A.size ();
182 dims(1) = 1;
183 Array<double> myb (dims);
184
185 for (std::size_t i = 0; i < A.size (); ++i)
186 myb.xelem (i) = A[i];
187
188 retval(0) = myb;
189 } 140 }
190 141
191 else if (a.rank () == 0) 142 else if (a.rank () == 0)
192 { 143 {
193 double b = dolfin::assemble (a); 144 double b = dolfin::assemble (a);
217 const dolfin::Form & a = frm.get_form (); 168 const dolfin::Form & a = frm.get_form ();
218 a.check (); 169 a.check ();
219 170
220 if (a.rank () == 1) 171 if (a.rank () == 1)
221 { 172 {
173 femfenics_factory factory;
174
222 dolfin::Vector A; 175 dolfin::Vector A;
223 dolfin::assemble (A, a); 176 dolfin::assemble (A, a);
224 177
225 dolfin::Vector x (myx.length ()); 178 dolfin::Vector x (myx.length ());
226 for (std::size_t i = 0; i < myx.length (); ++i) 179 for (std::size_t i = 0; i < myx.length (); ++i)
245 } 198 }
246 else 199 else
247 error ("assemble: unknown argument type"); 200 error ("assemble: unknown argument type");
248 } 201 }
249 202
250 dim_vector dims; 203 retval(0) = factory.vector (A);
251 dims.resize (2); 204 retval(1) = factory.vector (x);
252 dims(0) = A.size ();
253 dims(1) = 1;
254 Array<double> myb (dims), myc (dims);
255
256 for (std::size_t i = 0; i < A.size (); ++i)
257 {
258 myb.xelem (i) = A[i];
259 myc.xelem (i) = x[i];
260 }
261
262 retval(0) = myb;
263 retval(1) = myc;
264 } 205 }
265 206
266 else 207 else
267 error ("assemble: unknown size"); 208 error ("assemble: unknown size");
268 } 209 }