9 #include "assemble_matrix_impl.h"
10 #include "assemble_scalar_impl.h"
11 #include "assemble_vector_impl.h"
12 #include <Eigen/Dense>
13 #include <Eigen/Sparse>
41 return fem::impl::assemble_scalar(M);
54 fem::impl::assemble_vector(b, L);
77 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
78 const std::vector<std::shared_ptr<
const Form<T>>>& a,
79 const std::vector<std::vector<std::shared_ptr<
const DirichletBC<T>>>>& bcs1,
80 const std::vector<Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>>&
84 fem::impl::apply_lifting(b, a, bcs1, x0, scale);
97 const std::function<
int(std::int32_t,
const std::int32_t*, std::int32_t,
98 const std::int32_t*,
const T*)>& mat_add,
108 std::vector<bool> dof_marker0, dof_marker1;
110 = map0->block_size() * (map0->size_local() + map0->num_ghosts());
112 = map1->block_size() * (map1->size_local() + map1->num_ghosts());
113 for (std::size_t k = 0; k < bcs.size(); ++k)
116 assert(bcs[k]->function_space());
119 dof_marker0.resize(dim0,
false);
120 bcs[k]->mark_dofs(dof_marker0);
124 dof_marker1.resize(dim1,
false);
125 bcs[k]->mark_dofs(dof_marker1);
130 impl::assemble_matrix(mat_add, a, dof_marker0, dof_marker1);
143 template <
typename T>
145 const std::function<
int(std::int32_t,
const std::int32_t*, std::int32_t,
146 const std::int32_t*,
const T*)>& mat_add,
147 const Form<T>& a,
const std::vector<bool>& dof_marker0,
148 const std::vector<bool>& dof_marker1)
151 impl::assemble_matrix(mat_add, a, dof_marker0, dof_marker1);
160 template <
typename T>
166 std::vector<Eigen::Triplet<T>> triplets;
168 = [&triplets](std::int32_t nrow,
const std::int32_t* rows,
169 std::int32_t ncol,
const std::int32_t* cols,
const T* v) {
170 for (
int i = 0; i < nrow; ++i)
171 for (
int j = 0; j < ncol; ++j)
172 triplets.emplace_back(rows[i], cols[j], v[i * ncol + j]);
177 assemble_matrix<T>(mat_add, a, bcs);
181 Eigen::SparseMatrix<T, Eigen::RowMajor> mat(
182 map0->block_size() * (map0->size_local() + map0->num_ghosts()),
183 map1->block_size() * (map1->size_local() + map1->num_ghosts()));
184 mat.setFromTriplets(triplets.begin(), triplets.end());
198 template <
typename T>
200 const std::function<
int(std::int32_t,
const std::int32_t*, std::int32_t,
201 const std::int32_t*,
const T*)>& mat_add,
202 const Eigen::Ref<
const Eigen::Array<std::int32_t, Eigen::Dynamic, 1>>& rows,
205 for (Eigen::Index i = 0; i < rows.size(); ++i)
207 const std::int32_t row = rows(i);
208 mat_add(1, &row, 1, &row, &diagonal);
227 template <
typename T>
229 const std::function<
int(std::int32_t,
const std::int32_t*, std::int32_t,
230 const std::int32_t*,
const T*)>& mat_add,
235 for (
const auto& bc : bcs)
238 if (V.
contains(*bc->function_space()))
239 add_diagonal<T>(mat_add, bc->dofs_owned().col(0), diagonal);
253 template <
typename T>
254 void set_bc(Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
256 const Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>& x0,
259 if (b.rows() > x0.rows())
260 throw std::runtime_error(
"Size mismatch between b and x0 vectors.");
261 for (
const auto& bc : bcs)
264 bc->set(b, x0, scale);
271 template <
typename T>
272 void set_bc(Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
276 for (
const auto& bc : bcs)
293 template <
typename T>
294 std::vector<std::vector<std::shared_ptr<const fem::DirichletBC<T>>>>
299 std::vector<std::vector<std::shared_ptr<const fem::DirichletBC<T>>>> bcs0(
301 for (std::size_t i = 0; i < L.size(); ++i)
303 if (L[i]->function_space(0)->contains(*bc->function_space()))
304 bcs0[i].push_back(bc);
318 template <
typename T>
320 std::vector<std::vector<std::shared_ptr<const fem::DirichletBC<T>>>>>
326 std::vector<std::vector<std::shared_ptr<const fem::DirichletBC<T>>>>>
328 for (std::size_t i = 0; i < a.size(); ++i)
330 for (std::size_t j = 0; j < a[i].size(); ++j)
332 bcs1[i].resize(a[j].size());
338 if (a[i][j]->function_space(1)->contains(*bc->function_space()))
339 bcs1[i][j].push_back(bc);