27 : _dof_table_single_component(dof_table)
39 "The d.o.f. table passed must be for one variable that has only "
45 const unsigned num_components,
48 std::vector<GlobalVector*>
const& x,
49 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_table)
51 auto const num_nodal_dof_result =
54 std::vector<GlobalIndexType> ghost_indices;
58 auto const& single_component_ghost_indices =
60 auto const single_component_ghost_indices_size =
61 single_component_ghost_indices.size();
62 ghost_indices.reserve(single_component_ghost_indices_size *
64 for (
unsigned i = 0; i < single_component_ghost_indices_size; ++i)
66 for (
unsigned c = 0;
c < num_components; ++
c)
68 ghost_indices.push_back(
69 single_component_ghost_indices[i] * num_components +
c);
83 {num_nodal_dof_result, num_nodal_dof_result, &ghost_indices,
94 auto const size = extrapolatables.
size();
95 for (std::size_t i = 0; i < size; ++i)
107 const unsigned num_components,
110 std::vector<GlobalVector*>
const& x,
111 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_table)
125 if (
static_cast<std::size_t
>(num_element_dof_result) !=
126 extrapolatables.
size() * num_components)
128 OGS_FATAL(
"mismatch in number of D.o.F.");
131 auto const size = extrapolatables.
size();
132 for (std::size_t i = 0; i < size; ++i)
141 std::size_t
const element_index,
142 const unsigned num_components,
145 std::vector<GlobalVector*>
const& x,
146 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_table,
149 auto const& integration_point_values =
154 if (integration_point_values.empty())
159 auto const& N_0 = extrapolatables.
getShapeMatrix(element_index, 0);
160 auto const num_nodes =
static_cast<unsigned>(N_0.cols());
161 auto const num_values =
162 static_cast<unsigned>(integration_point_values.size());
164 if (num_values % num_components != 0)
167 "The number of computed integration point values is not divisible "
168 "by the number of num_components. Maybe the computed property is "
169 "not a {:d}-component vector for each integration point.",
174 const auto num_int_pts = num_values / num_components;
176 if (num_int_pts < num_nodes)
179 "Least squares is not possible if there are more nodes than "
180 "integration points.");
184 std::make_pair(num_nodes, num_int_pts),
CachedData{});
186 auto& cached_data = pair_it_inserted.first->second;
187 if (pair_it_inserted.second)
189 DBUG(
"Computing new singular value decomposition");
194 auto& interpolation_matrix = cached_data.A;
195 interpolation_matrix.resize(num_int_pts, num_nodes);
197 interpolation_matrix.row(0) = N_0;
198 for (
unsigned int_pt = 1; int_pt < num_int_pts; ++int_pt)
200 auto const& shp_mat =
202 assert(shp_mat.cols() == num_nodes);
205 interpolation_matrix.row(int_pt) = shp_mat;
214 Eigen::JacobiSVD<Eigen::MatrixXd> svd(
215 interpolation_matrix, Eigen::ComputeThinU | Eigen::ComputeThinV);
217 auto const& S = svd.singularValues();
218 auto const& U = svd.matrixU();
219 auto const& V = svd.matrixV();
222 auto const rank = svd.rank();
223 assert(rank == num_nodes);
226 cached_data.A_pinv.noalias() = V.leftCols(rank) *
227 S.head(rank).asDiagonal().inverse() *
228 U.leftCols(rank).transpose();
230 else if (cached_data.A.row(0) != N_0)
232 OGS_FATAL(
"The cached and the passed shapematrices differ.");
235 auto const& global_indices =
238 if (num_components == 1)
240 auto const integration_point_values_vec =
244 Eigen::VectorXd
const nodal_values =
245 cached_data.A_pinv * integration_point_values_vec;
250 counts.
add(global_indices,
251 std::vector<double>(global_indices.size(), 1.0));
256 integration_point_values, num_components, num_int_pts);
259 Eigen::MatrixXd
const nodal_values =
260 cached_data.A_pinv * integration_point_values_mat.transpose();
262 std::vector<GlobalIndexType> indices;
263 indices.reserve(num_components * global_indices.size());
266 for (
unsigned comp = 0; comp < num_components; ++comp)
268 transform(cbegin(global_indices), cend(global_indices),
269 back_inserter(indices),
270 [&](
auto const i) {
return num_components * i + comp; });
276 counts.
add(indices, std::vector<double>(indices.size(), 1.0));
281 std::size_t
const element_index,
282 const unsigned num_components,
285 std::vector<GlobalVector*>
const& x,
286 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_table)
291 auto const num_values =
static_cast<unsigned>(int_pt_vals.size());
292 if (num_values % num_components != 0)
295 "The number of computed integration point values is not divisible "
296 "by the number of num_components. Maybe the computed property is "
297 "not a {:d}-component vector for each integration point.",
302 const auto num_int_pts = num_values / num_components;
304 const auto& global_indices =
306 const auto num_nodes =
static_cast<unsigned>(global_indices.size());
308 auto const& interpolation_matrix =
311 Eigen::VectorXd nodal_vals_element(num_nodes);
312 auto const int_pt_vals_mat =
317 for (
unsigned comp = 0; comp < num_components; ++comp)
320 for (
unsigned i = 0; i < num_nodes; ++i)
323 auto const idx = num_components * global_indices[i] + comp;
327 double const residual = (interpolation_matrix * nodal_vals_element -
328 int_pt_vals_mat.row(comp).transpose())
334 auto const root_mean_square = std::sqrt(residual / num_int_pts);
GlobalMatrix::IndexType GlobalIndexType
void DBUG(char const *fmt, Args const &... args)
Global vector based on Eigen vector.
void add(IndexType rowId, double v)
add entry
std::vector< GlobalIndexType > const & getGhostIndices() const
Get ghost indices, forwarded from MeshComponentMap.
std::size_t dofSizeWithoutGhosts() const
int getNumberOfGlobalComponents() const
void finalizeAssembly(PETScMatrix &A)
void componentwiseDivide(PETScVector &w, PETScVector const &x, PETScVector const &y)
void setLocalAccessibleVector(PETScVector const &x)
Eigen::Map< const Vector > toVector(std::vector< double > const &data, Eigen::VectorXd::Index size)
Creates an Eigen mapped vector from the given data vector.
Eigen::Map< const Matrix > toMatrix(std::vector< double > const &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)