1 #include <Xpetra_MatrixFactory.hpp>
2 #include <Xpetra_Matrix.hpp>
3 #include <Xpetra_MatrixMatrix.hpp>
4 #include <Xpetra_CrsMatrixWrap.hpp>
5 #include <Xpetra_BlockedCrsMatrix.hpp>
6 #include <Xpetra_CrsMatrix.hpp>
7 #include <MueLu_Utilities_decl.hpp>
27 template <
typename scalar>
29 Teuchos::RCP<Teuchos::ParameterList> randomParams,
30 Kokkos::View<scalar *, Kokkos::LayoutLeft> m_rv_values)
32 numPrimalDPN = _numPrimalDPN;
35 Teuchos::RCP<Map> map = Teuchos::rcp(
new Map(pbl, numPrimalDPN));
38 domain = Teuchos::rcp(
new Domain<scalar>(pbl, map, randomParams, m_rv_values));
55 template <
typename scalar>
62 const size_t myRank = algebraic->map->mapDofs->getComm()->getRank();
66 RCP<tpetra_mvector_type> mx =
67 (Teuchos::rcp_dynamic_cast<xpetra_tmvector_type>(algebraic->vectors->solutionMultiVector))->getTpetra_MultiVector();
68 RCP<tpetra_vector_type> xl = mx->getVectorNonConst(0);
70 RCP<tpetra_vector_type> x = rcp(
new tpetra_vector_type(algebraic->map->mapDofs,
false));
72 RCP<xpetra_bmvector_type> blockedSol =
73 rcp(
new xpetra_bmvector_type(algebraic->map->blockedMap, algebraic->vectors->solutionMultiVector));
75 RCP<tpetra_vector_type> l =
76 (Teuchos::rcp_dynamic_cast<xpetra_tmvector_type>(blockedSol->getMultiVector(1,
true)))->getTpetra_MultiVector()->getVectorNonConst(0);
78 RCP<tpetra_vector_type> gap = rcp(
new tpetra_vector_type(*algebraic->vectors->initialGap, Teuchos::Copy));
80 RCP<export_type> exportx = rcp(
new export_type(algebraic->map->fullmap, algebraic->map->mapDofs));
82 x->doExport(*xl, *exportx, Tpetra::INSERT);
88 algebraic->matrices->B_G->apply(*x, *gap, Teuchos::NO_TRANS, ((scalar)1.), ((scalar)-1.));
92 l->update((scalar)c, *gap, (scalar)1.);
96 l->template sync<Kokkos::HostSpace>();
97 auto l_2d = l->template getLocalView<Kokkos::HostSpace>();
99 Mask<scalar> hasConverged;
100 const size_t ensemble_size = hasConverged.getSize();
102 for (
size_t l = 0; l < ensemble_size; ++l)
103 hasConverged.set(l,
true);
105 bool hasConverged_old =
false;
106 if (active_set_iteration > 1)
107 hasConverged_old =
true;
109 algebraic->matrices->B_2->resumeFill();
110 algebraic->matrices->C->resumeFill();
111 algebraic->matrices->Cb->resumeFill();
113 const size_t numDPN = algebraic->map->numPrimalDPN;
116 if (domain->contactsList->hasAtLeastOneSticking() && numDPN == 4)
118 else if (domain->contactsList->hasAtLeastOneSticking())
120 else if (numDPN == 4)
125 const int numMyLagrangeNum = algebraic->map->mapLagrangeDofs->getNodeNumElements();
127 const int numMyMechLagrangeNum = numMyLagrangeNum / numLMPN;
129 for (
auto local_index = 0; local_index < numMyMechLagrangeNum; ++local_index)
131 global_ordinal_type lagrange_id = algebraic->map->mapLagrangeDofs->getGlobalElement(local_index);
133 lagrange_id = lagrange_id * numLMPN;
134 local_ordinal_type interface_id = domain->contactsList->getInterfaceOfSlaveNode(node_id);
136 if (domain->contactsList->isNotUpdated(interface_id))
139 scalar old_mask = old_activity(local_index);
140 scalar old_old_mask = old_old_activity(local_index);
141 Mask<scalar> current_mask = (l_2d(local_index * numLMPN, 0) >= 0.);
143 scalar oneMask, oneNotMask;
144 mask_assign<scalar>(current_mask, oneMask) = {1., 0.};
145 mask_assign<scalar>(current_mask, oneNotMask) = {0., 1.};
147 if (active_set_iteration > 1)
148 if (!MaskLogic::AND(oneMask == old_old_mask))
149 hasConverged_old =
false;
151 hasConverged = hasConverged && (oneMask == old_mask);
153 if (!MaskLogic::AND(oneMask == old_mask))
155 global_ordinal_type lagrange_id = algebraic->map->mapLagrangeDofs->getGlobalElement(local_index);
157 lagrange_id = lagrange_id * numLMPN;
161 if (MaskLogic::AND(current_mask))
163 Teuchos::ArrayView<const local_ordinal_type> local_indices;
164 Teuchos::ArrayView<const scalar> values;
166 if (domain->contactsList->isSticking(interface_id))
168 for (
size_t j = 0; j < numLMPN; ++j)
170 algebraic->matrices->B->getLocalRowView(local_index * numLMPN + j, local_indices, values);
172 for (
auto k = 0; k < local_indices.size(); ++k)
174 global_ordinal_type column_index = algebraic->matrices->B->getColMap()->getGlobalElement(local_indices[k]);
176 algebraic->matrices->B_2->replaceGlobalValues(lagrange_id + j,
177 Teuchos::tuple<global_ordinal_type>(column_index),
178 Teuchos::tuple<scalar>(values[k]));
181 algebraic->matrices->C->replaceGlobalValues(lagrange_id + j,
182 Teuchos::tuple<global_ordinal_type>(lagrange_id + j),
183 Teuchos::tuple<scalar>((scalar)0.));
185 algebraic->matrices->Cb->replaceGlobalValues(lagrange_id + j,
186 Teuchos::tuple<global_ordinal_type>(lagrange_id + j),
187 Teuchos::tuple<scalar>((scalar)1.));
192 algebraic->matrices->B->getLocalRowView(local_index * numLMPN, local_indices, values);
194 for (
auto k = 0; k < local_indices.size(); ++k)
196 global_ordinal_type column_index = algebraic->matrices->B->getColMap()->getGlobalElement(local_indices[k]);
198 algebraic->matrices->B_2->replaceGlobalValues(lagrange_id,
199 Teuchos::tuple<global_ordinal_type>(column_index),
200 Teuchos::tuple<scalar>(values[k]));
203 algebraic->matrices->C->replaceGlobalValues(lagrange_id,
204 Teuchos::tuple<global_ordinal_type>(lagrange_id),
205 Teuchos::tuple<scalar>((scalar)0.));
207 algebraic->matrices->Cb->replaceGlobalValues(lagrange_id,
208 Teuchos::tuple<global_ordinal_type>(lagrange_id),
209 Teuchos::tuple<scalar>((scalar)1.));
211 if (numLMPN == 2 || numLMPN == 4)
213 algebraic->matrices->B->getLocalRowView(local_index * numLMPN + numLMPN - 1, local_indices, values);
215 for (
auto k = 0; k < local_indices.size(); ++k)
217 global_ordinal_type column_index = algebraic->matrices->B->getColMap()->getGlobalElement(local_indices[k]);
219 algebraic->matrices->B_2->replaceGlobalValues(lagrange_id + numLMPN - 1,
220 Teuchos::tuple<global_ordinal_type>(column_index),
221 Teuchos::tuple<scalar>(values[k]));
224 algebraic->matrices->C->replaceGlobalValues(lagrange_id + numLMPN - 1,
225 Teuchos::tuple<global_ordinal_type>(lagrange_id + numLMPN - 1),
226 Teuchos::tuple<scalar>((scalar)0.));
228 algebraic->matrices->Cb->replaceGlobalValues(lagrange_id + numLMPN - 1,
229 Teuchos::tuple<global_ordinal_type>(lagrange_id + numLMPN - 1),
230 Teuchos::tuple<scalar>((scalar)1.));
236 else if (!MaskLogic::AND(current_mask))
238 Teuchos::ArrayView<const local_ordinal_type> local_indices;
239 Teuchos::ArrayView<const scalar> values;
241 if (domain->contactsList->isSticking(interface_id))
243 for (
size_t j = 0; j < numLMPN; ++j)
245 algebraic->matrices->B->getLocalRowView(local_index * numLMPN + j, local_indices, values);
247 for (
auto k = 0; k < local_indices.size(); ++k)
249 global_ordinal_type column_index = algebraic->matrices->B->getColMap()->getGlobalElement(local_indices[k]);
251 algebraic->matrices->B_2->replaceGlobalValues(lagrange_id + j,
252 Teuchos::tuple<global_ordinal_type>(column_index),
253 Teuchos::tuple<scalar>((scalar)0.));
256 algebraic->matrices->C->replaceGlobalValues(lagrange_id + j,
257 Teuchos::tuple<global_ordinal_type>(lagrange_id + j),
258 Teuchos::tuple<scalar>((scalar)1.));
260 algebraic->matrices->Cb->replaceGlobalValues(lagrange_id + j,
261 Teuchos::tuple<global_ordinal_type>(lagrange_id + j),
262 Teuchos::tuple<scalar>((scalar)0.));
267 algebraic->matrices->B->getLocalRowView(local_index * numLMPN, local_indices, values);
269 for (
auto k = 0; k < local_indices.size(); ++k)
271 global_ordinal_type column_index = algebraic->matrices->B->getColMap()->getGlobalElement(local_indices[k]);
273 algebraic->matrices->B_2->replaceGlobalValues(lagrange_id,
274 Teuchos::tuple<global_ordinal_type>(column_index),
275 Teuchos::tuple<scalar>((scalar)0.));
278 algebraic->matrices->C->replaceGlobalValues(lagrange_id,
279 Teuchos::tuple<global_ordinal_type>(lagrange_id),
280 Teuchos::tuple<scalar>((scalar)1.));
282 algebraic->matrices->Cb->replaceGlobalValues(lagrange_id,
283 Teuchos::tuple<global_ordinal_type>(lagrange_id),
284 Teuchos::tuple<scalar>((scalar)0.));
286 if (numLMPN == 2 || numLMPN == 4)
288 algebraic->matrices->B->getLocalRowView(local_index * numLMPN + numLMPN - 1, local_indices, values);
290 for (
auto k = 0; k < local_indices.size(); ++k)
292 global_ordinal_type column_index = algebraic->matrices->B->getColMap()->getGlobalElement(local_indices[k]);
294 algebraic->matrices->B_2->replaceGlobalValues(lagrange_id + numLMPN - 1,
295 Teuchos::tuple<global_ordinal_type>(column_index),
296 Teuchos::tuple<scalar>((scalar)0.));
299 algebraic->matrices->C->replaceGlobalValues(lagrange_id + numLMPN - 1,
300 Teuchos::tuple<global_ordinal_type>(lagrange_id + numLMPN - 1),
301 Teuchos::tuple<scalar>((scalar)1.));
303 algebraic->matrices->Cb->replaceGlobalValues(lagrange_id + numLMPN - 1,
304 Teuchos::tuple<global_ordinal_type>(lagrange_id + numLMPN - 1),
305 Teuchos::tuple<scalar>((scalar)0.));
313 Teuchos::ArrayView<const local_ordinal_type> local_indices;
314 Teuchos::ArrayView<const scalar> values;
316 if (domain->contactsList->isSticking(interface_id))
318 for (
size_t j = 0; j < numLMPN; ++j)
320 algebraic->matrices->B->getLocalRowView(local_index * numLMPN + j, local_indices, values);
322 for (
auto k = 0; k < local_indices.size(); ++k)
324 global_ordinal_type column_index = algebraic->matrices->B->getColMap()->getGlobalElement(local_indices[k]);
326 algebraic->matrices->B_2->replaceGlobalValues(lagrange_id + j,
327 Teuchos::tuple<global_ordinal_type>(column_index),
328 Teuchos::tuple<scalar>(oneMask * values[k]));
331 algebraic->matrices->C->replaceGlobalValues(lagrange_id + j,
332 Teuchos::tuple<global_ordinal_type>(lagrange_id + j),
333 Teuchos::tuple<scalar>(oneNotMask));
335 algebraic->matrices->Cb->replaceGlobalValues(lagrange_id + j,
336 Teuchos::tuple<global_ordinal_type>(lagrange_id + j),
337 Teuchos::tuple<scalar>(oneMask));
342 algebraic->matrices->B->getLocalRowView(local_index * numLMPN, local_indices, values);
344 for (
auto k = 0; k < local_indices.size(); ++k)
346 global_ordinal_type column_index = algebraic->matrices->B->getColMap()->getGlobalElement(local_indices[k]);
348 algebraic->matrices->B_2->replaceGlobalValues(lagrange_id,
349 Teuchos::tuple<global_ordinal_type>(column_index),
350 Teuchos::tuple<scalar>(oneMask * values[k]));
353 algebraic->matrices->C->replaceGlobalValues(lagrange_id,
354 Teuchos::tuple<global_ordinal_type>(lagrange_id),
355 Teuchos::tuple<scalar>(oneNotMask));
357 algebraic->matrices->Cb->replaceGlobalValues(lagrange_id,
358 Teuchos::tuple<global_ordinal_type>(lagrange_id),
359 Teuchos::tuple<scalar>(oneMask));
361 if (numLMPN == 2 || numLMPN == 4)
363 algebraic->matrices->B->getLocalRowView(local_index * numLMPN + numLMPN - 1, local_indices, values);
365 for (
auto k = 0; k < local_indices.size(); ++k)
367 global_ordinal_type column_index = algebraic->matrices->B->getColMap()->getGlobalElement(local_indices[k]);
369 algebraic->matrices->B_2->replaceGlobalValues(lagrange_id + numLMPN - 1,
370 Teuchos::tuple<global_ordinal_type>(column_index),
371 Teuchos::tuple<scalar>(oneMask * values[k]));
374 algebraic->matrices->C->replaceGlobalValues(lagrange_id + numLMPN - 1,
375 Teuchos::tuple<global_ordinal_type>(lagrange_id + numLMPN - 1),
376 Teuchos::tuple<scalar>(oneNotMask));
378 algebraic->matrices->Cb->replaceGlobalValues(lagrange_id + numLMPN - 1,
379 Teuchos::tuple<global_ordinal_type>(lagrange_id + numLMPN - 1),
380 Teuchos::tuple<scalar>(oneMask));
386 old_old_activity(local_index) = old_activity(local_index);
387 old_activity(local_index) = oneMask;
390 algebraic->matrices->B_2->fillComplete(algebraic->map->mapDofs, algebraic->map->mapLagrangeDofs);
391 algebraic->matrices->C->fillComplete();
392 algebraic->matrices->Cb->fillComplete();
394 RCP<tpetra_vector_type> initialActiveGap(
new tpetra_vector_type(algebraic->map->mapLagrangeDofs,
true));
395 algebraic->matrices->Cb->apply(*algebraic->vectors->initialGap, *initialActiveGap);
397 RCP<xpetra_mvector_type> xg =
399 algebraic->vectors->rhsBlockedMultiVector->setMultiVector(1, xg,
true);
400 algebraic->vectors->rhsMultiVector = algebraic->vectors->rhsBlockedMultiVector->Merge();
402 size_t numtasks = algebraic->map->mapDofs->getComm()->getSize();
408 Mask<scalar> globalHasConverged = hasConverged;
409 bool globalHasConverged_old = hasConverged_old;
411 for (
size_t i = 1; i < numtasks; ++i)
413 for (
size_t l = 0; l < ensemble_size; ++l)
415 MPI_Recv(&tmp, 1, MPI_DOUBLE, i, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
417 globalHasConverged.set(l,
false);
419 MPI_Recv(&tmp, 1, MPI_DOUBLE, i, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
421 globalHasConverged_old =
false;
423 hasConverged = globalHasConverged;
424 hasConverged_old = globalHasConverged_old;
426 for (
size_t i = 1; i < numtasks; ++i)
428 for (
size_t l = 0; l < ensemble_size; ++l)
430 tmp = hasConverged.get(l) ? 1. : -1.;
431 MPI_Send(&tmp, 1, MPI_DOUBLE, i, 0, MPI_COMM_WORLD);
433 tmp = hasConverged_old ? 1. : -1.;
434 MPI_Send(&tmp, 1, MPI_DOUBLE, i, 0, MPI_COMM_WORLD);
439 for (
size_t l = 0; l < ensemble_size; ++l)
441 tmp = hasConverged.get(l) ? 1. : -1.;
442 MPI_Send(&tmp, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
444 tmp = hasConverged_old ? 1. : -1.;
445 MPI_Send(&tmp, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
446 for (
size_t l = 0; l < ensemble_size; ++l)
448 MPI_Recv(&tmp, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
449 hasConverged.set(l, tmp > 0.);
451 MPI_Recv(&tmp, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
452 hasConverged_old = (tmp > 0.);
455 if (hasConverged_old)
459 std::cout <<
"-----------------------" << std::endl;
460 std::cout <<
"---Has not converged---" << std::endl;
461 std::cout <<
"-----------------------" << std::endl;
478 template <
typename scalar>
481 using Teuchos::Array;
485 const std::vector<RCP<const tpetra_map_type>> maps = {algebraic->map->mapDofs, algebraic->map->mapLagrangeDofs};
487 const size_t numDofs = algebraic->map->mapDofs->getGlobalNumElements();
488 const size_t numMyDofs = algebraic->map->mapDofs->getNodeNumElements();
489 const size_t numMyDofsWO = algebraic->map->mapDofsWO->getNodeNumElements();
490 Teuchos::ArrayView<const global_ordinal_type> myDofs = algebraic->map->mapDofs->getNodeElementList();
491 Teuchos::ArrayView<const global_ordinal_type> myDofsWO = algebraic->map->mapDofsWO->getNodeElementList();
493 const size_t numLagrangeDofs = algebraic->map->mapLagrangeDofs->getGlobalNumElements();
494 const size_t numMyLagrangeDofs = algebraic->map->mapLagrangeDofs->getNodeNumElements();
495 Teuchos::ArrayView<const global_ordinal_type> myLagrangeDofs = algebraic->map->mapLagrangeDofs->getNodeElementList();
497 Array<global_ordinal_type> myFullDofs(numMyDofs + numMyLagrangeDofs);
498 for (
auto i = 0; i < numMyDofs; ++i)
499 myFullDofs[i] = myDofs[i];
500 for (
auto i = 0; i < numMyLagrangeDofs; ++i)
501 myFullDofs[numMyDofs + i] = numDofs + myLagrangeDofs[i];
503 Array<global_ordinal_type> myFullDofsWO(numMyDofsWO + numMyLagrangeDofs);
504 for (
auto i = 0; i < numMyDofsWO; ++i)
505 myFullDofsWO[i] = myDofsWO[i];
506 for (
auto i = 0; i < numMyLagrangeDofs; ++i)
507 myFullDofsWO[numMyDofsWO + i] = numDofs + myLagrangeDofs[i];
509 algebraic->map->fullmap =
512 algebraic->map->mapDofs->getIndexBase(),
513 algebraic->map->mapDofs->getComm()));
514 algebraic->map->fullmapWO =
517 algebraic->map->mapDofs->getIndexBase(),
518 algebraic->map->mapDofs->getComm()));
520 RCP<xpetra_map_type> xfullmap = rcp(
new xpetra_tmap_type(algebraic->map->fullmap));
522 const size_t myRank = algebraic->map->mapDofs->getComm()->getRank();
523 const size_t numLclIndices = (myRank == 0) ? numDofs + numLagrangeDofs : 0;
525 algebraic->map->fullmapOutput =
528 algebraic->map->mapDofs->getIndexBase(),
529 algebraic->map->mapDofs->getComm()));
531 std::vector<RCP<const xpetra_map_type>> xsubmaps = {rcp(
new xpetra_tmap_type(algebraic->map->mapDofs)),
534 algebraic->map->blockedMap =
537 algebraic->matrices->blockedMatrix =
538 rcp(
new xpetra_bcrs_type(algebraic->map->blockedMap, algebraic->map->blockedMap, 8));
542 algebraic->matrices->B_2->resumeFill();
543 algebraic->matrices->C->resumeFill();
544 algebraic->matrices->Cb->resumeFill();
546 const size_t numDPN = algebraic->map->numPrimalDPN;
549 if (domain->contactsList->hasAtLeastOneSticking() && numDPN == 4)
551 else if (domain->contactsList->hasAtLeastOneSticking())
553 else if (numDPN == 4)
558 const int numMyLagrangeNum = algebraic->map->mapLagrangeDofs->getNodeNumElements();
560 const int numMyMechLagrangeNum = numMyLagrangeNum / numLMPN;
562 old_activity = Kokkos::View<scalar *, Kokkos::LayoutRight>(
"R", numMyMechLagrangeNum);
563 old_old_activity = Kokkos::View<scalar *, Kokkos::LayoutRight>(
"R", numMyMechLagrangeNum);
565 for (
auto local_index = 0; local_index < numMyMechLagrangeNum; ++local_index)
567 global_ordinal_type lagrange_id = algebraic->map->mapLagrangeDofs->getGlobalElement(local_index);
569 lagrange_id = lagrange_id * numLMPN;
570 local_ordinal_type interface_id = domain->contactsList->getInterfaceOfSlaveNode(node_id);
571 if (domain->contactsList->isInitiallyOpen(interface_id))
573 if (domain->contactsList->isNodeInitiallyClosed(interface_id, node_id))
575 old_activity(local_index) = (scalar)1.;
578 old_activity(local_index) = (scalar)0.;
580 Teuchos::ArrayView<const local_ordinal_type> local_indices;
581 Teuchos::ArrayView<const scalar> values;
583 for (
size_t j = 0; j < numLMPN; ++j)
585 algebraic->matrices->B_2->getLocalRowView(local_index * numLMPN + j, local_indices, values);
587 for (
auto k = 0; k < local_indices.size(); ++k)
589 global_ordinal_type column_index = algebraic->matrices->B->getColMap()->getGlobalElement(local_indices[k]);
591 algebraic->matrices->B_2->replaceGlobalValues(lagrange_id + j,
592 Teuchos::tuple<global_ordinal_type>(column_index),
593 Teuchos::tuple<scalar>((scalar)0.));
596 algebraic->matrices->C->replaceGlobalValues(lagrange_id + j,
597 Teuchos::tuple<global_ordinal_type>(lagrange_id + j),
598 Teuchos::tuple<scalar>((scalar)1.));
600 algebraic->matrices->Cb->replaceGlobalValues(lagrange_id + j,
601 Teuchos::tuple<global_ordinal_type>(lagrange_id + j),
602 Teuchos::tuple<scalar>((scalar)0.));
607 old_activity(local_index) = (scalar)1.;
608 if (domain->contactsList->isSticking(interface_id) ==
false && domain->contactsList->hasAtLeastOneSticking() ==
true)
610 Teuchos::ArrayView<const local_ordinal_type> local_indices;
611 Teuchos::ArrayView<const scalar> values;
613 for (
size_t j = 1; j < numLMPN; ++j)
615 algebraic->matrices->B_2->getLocalRowView(local_index * numLMPN + j, local_indices, values);
617 for (
auto k = 0; k < local_indices.size(); ++k)
619 global_ordinal_type column_index = algebraic->matrices->B->getColMap()->getGlobalElement(local_indices[k]);
621 algebraic->matrices->B_2->replaceGlobalValues(lagrange_id + j,
622 Teuchos::tuple<global_ordinal_type>(column_index),
623 Teuchos::tuple<scalar>((scalar)0.));
626 algebraic->matrices->C->replaceGlobalValues(lagrange_id + j,
627 Teuchos::tuple<global_ordinal_type>(lagrange_id + j),
628 Teuchos::tuple<scalar>((scalar)1.));
630 algebraic->matrices->Cb->replaceGlobalValues(lagrange_id + j,
631 Teuchos::tuple<global_ordinal_type>(lagrange_id + j),
632 Teuchos::tuple<scalar>((scalar)0.));
637 algebraic->matrices->B_2->fillComplete(algebraic->map->mapDofs, algebraic->map->mapLagrangeDofs);
638 algebraic->matrices->C->fillComplete();
639 algebraic->matrices->Cb->fillComplete();
641 algebraic->matrices->xB = rcp(
new xpetra_tcrs_type(algebraic->matrices->B_2));
642 algebraic->matrices->xB_T = rcp(
new xpetra_tcrs_type(algebraic->matrices->B_T));
643 algebraic->matrices->xC = rcp(
new xpetra_tcrs_type(algebraic->matrices->C));
645 algebraic->matrices->xwB = rcp(
new xpetra_wcrs_type(algebraic->matrices->xB));
646 algebraic->matrices->xwB_T = rcp(
new xpetra_wcrs_type(algebraic->matrices->xB_T));
647 algebraic->matrices->xwC = rcp(
new xpetra_wcrs_type(algebraic->matrices->xC));
649 RCP<Xpetra::Matrix<scalar, local_ordinal_type, global_ordinal_type, node_type>> CTC = Teuchos::null;
650 RCP<Xpetra::Matrix<scalar, local_ordinal_type, global_ordinal_type, node_type>> ApCTC = Teuchos::null;
651 std::string MatrixMatrix_output_file =
"MatrixMatrix_out.txt";
652 std::ofstream MatrixMatrix_ofstream(MatrixMatrix_output_file);
653 RCP<Teuchos::FancyOStream> verbOut = Teuchos::getFancyOStream(Teuchos::rcpFromRef(MatrixMatrix_ofstream));
655 Xpetra::MatrixMatrix<scalar, local_ordinal_type, global_ordinal_type, node_type>::Multiply(*(algebraic->matrices->xwB_T),
657 *(algebraic->matrices->xwB),
660 RCP<xpetra_bcrs_type> bA = Teuchos::rcp_dynamic_cast<xpetra_bcrs_type>(algebraic->matrices->A);
664 if (bA != Teuchos::null && merge)
665 Xpetra::MatrixMatrix<scalar, local_ordinal_type, global_ordinal_type, node_type>::TwoMatrixAdd(*CTC,
674 Xpetra::MatrixMatrix<scalar, local_ordinal_type, global_ordinal_type, node_type>::TwoMatrixAdd(*CTC,
677 *(algebraic->matrices->A),
685 if (bA != Teuchos::null && merge)
688 ApCTC = algebraic->matrices->A;
691 algebraic->matrices->blockedMatrix->setMatrix(0, 0, ApCTC);
692 algebraic->matrices->blockedMatrix->setMatrix(0, 1, algebraic->matrices->xwB_T);
693 algebraic->matrices->blockedMatrix->setMatrix(1, 0, algebraic->matrices->xwB);
694 algebraic->matrices->blockedMatrix->setMatrix(1, 1, algebraic->matrices->xwC);
696 algebraic->matrices->blockedMatrix->fillComplete();
697 RCP<const xpetra_wcrs_type> tmpOp = Teuchos::rcp_dynamic_cast<const xpetra_wcrs_type>(ApCTC);
698 const RCP<const xpetra_tcrs_type> &tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const xpetra_tcrs_type>(tmpOp->getCrsMatrix());
699 RCP<const tpetra_crs_type> tpetra_ApCTC = tmp_ECrsMtx->getTpetra_CrsMatrix();
702 Tpetra::MatrixMarket::Writer<tpetra_crs_type>::writeSparseFile(std::string(
"ApCTC_mm.txt"), tpetra_ApCTC);
713 template <
typename scalar>
716 using Teuchos::Array;
720 RCP<tpetra_vector_type> intialContactPosition(
new tpetra_vector_type(algebraic->map->mapDofs,
true));
721 algebraic->vectors->initialGap = rcp(
new tpetra_vector_type(algebraic->map->mapLagrangeDofs,
true));
722 RCP<tpetra_vector_type> initialActiveGap(
new tpetra_vector_type(algebraic->map->mapLagrangeDofs,
true));
724 algebraic->vectors->lagrange = rcp(
new tpetra_vector_type(algebraic->map->mapLagrangeDofs,
true));
727 size_t numMyContacts = domain->contactsList->getContactNumber();
728 size_t numPrimalDPN = algebraic->map->numPrimalDPN;
729 size_t numMPrimalDPN = 3;
731 if (numPrimalDPN == 4)
732 offSet = algebraic->map->mapNodes->getGlobalNumElements();
734 for (
auto i = 0; i < numMyContacts; ++i)
736 size_t numMySlaveNodes = domain->contactsList->getSlaveNodesSize(i);
737 for (
auto j = 0; j < numMySlaveNodes; ++j)
740 auto global_nodeid = domain->contactsList->getSlaveNode(i, j);
741 if (algebraic->map->mapNodes->isNodeGlobalElement(global_nodeid))
743 auto local_nodeid = algebraic->map->mapNodesWO->getLocalElement(global_nodeid);
744 for (
auto k = 0; k < 3; ++k)
746 scalar tmp = domain->nodesList->nodes(local_nodeid, k) +
747 domain->contactsList->getSlaveNormal(i, k) *
748 loads->preloadList->getPreloadValue(i);
749 intialContactPosition->replaceGlobalValue(offSet + global_nodeid * numMPrimalDPN + k, tmp);
754 size_t numMyMasterNodes = domain->contactsList->getMasterNodesSize(i);
755 for (
auto j = 0; j < numMyMasterNodes; ++j)
758 auto global_nodeid = domain->contactsList->getMasterNode(i, j);
759 if (algebraic->map->mapNodes->isNodeGlobalElement(global_nodeid))
761 auto local_nodeid = algebraic->map->mapNodesWO->getLocalElement(global_nodeid);
762 for (
auto k = 0; k < 3; ++k)
763 intialContactPosition->replaceGlobalValue(offSet + global_nodeid * numMPrimalDPN + k,
764 ((scalar)domain->nodesList->nodes(local_nodeid, k)));
770 algebraic->matrices->B_G->apply(*intialContactPosition, *algebraic->vectors->initialGap, Teuchos::NO_TRANS, ((scalar)-1.));
772 algebraic->matrices->Cb->apply(*algebraic->vectors->initialGap, *initialActiveGap);
775 if (domain->contactsList->isTying(0))
776 initialActiveGap->scale((scalar)0.);
779 Tpetra::MatrixMarket::Writer<tpetra_crs_type>::writeDenseFile(std::string(
"initialActiveGap_mm.txt"), initialActiveGap);
781 RCP<xpetra_mvector_type> xb =
784 RCP<xpetra_mvector_type> xg =
787 RCP<xpetra_mvector_type> xx =
790 RCP<xpetra_mvector_type> xlagrange =
793 algebraic->vectors->rhsBlockedMultiVector =
796 algebraic->vectors->solutionBlockedMultiVector =
799 algebraic->vectors->rhsBlockedMultiVector->setMultiVector(0, xb,
true);
800 algebraic->vectors->rhsBlockedMultiVector->setMultiVector(1, xg,
true);
802 algebraic->vectors->solutionBlockedMultiVector->setMultiVector(0, xx,
true);
803 algebraic->vectors->solutionBlockedMultiVector->setMultiVector(1, xlagrange,
true);
805 algebraic->vectors->rhsMultiVector = algebraic->vectors->rhsBlockedMultiVector->Merge();
806 algebraic->vectors->solutionMultiVector = algebraic->vectors->solutionBlockedMultiVector->Merge();
817 template <
typename scalar>
820 compute_B<scalar>(algebraic->matrices,
821 domain->elementsList,
823 domain->contactsList,
824 algebraic->map, comm);
830 template <
typename scalar>
836 typedef KokkosSparse::CrsMatrix<scalar, local_ordinal_type, typename tpetra_crs_type::execution_space> local_matrix_type;
838 typedef Kokkos::View<scalar ***, Kokkos::LayoutRight, Kokkos::DefaultExecutionSpace> elem_matrix_view_type;
840 local_matrix_type local_L, local_K, local_S;
842 if (numPrimalDPN == 4 || numPrimalDPN == 1)
844 algebraic->matrices->L->resumeFill();
845 algebraic->matrices->L->setAllToScalar((scalar)0.0);
847 local_L = algebraic->matrices->L->getLocalMatrix();
849 if (numPrimalDPN == 4 || numPrimalDPN == 3)
851 algebraic->matrices->K->resumeFill();
852 algebraic->matrices->K->setAllToScalar((scalar)0.0);
854 local_K = algebraic->matrices->K->getLocalMatrix();
856 if (numPrimalDPN == 4)
858 algebraic->matrices->S->resumeFill();
859 algebraic->matrices->S->setAllToScalar((scalar)0.0);
861 local_S = algebraic->matrices->S->getLocalMatrix();
864 auto numMyElems = domain->elementsList->getElementNumber();
866 elem_matrix_view_type K, L, S;
868 #ifdef KOKKOS_ENABLE_DEPRECATED_CODE
869 int pool_size = Kokkos::DefaultExecutionSpace::thread_pool_size();
871 int pool_size = Kokkos::DefaultExecutionSpace::impl_thread_pool_size();
873 int max_element_size = 1;
875 for (
int i_elem = 0; i_elem < numMyElems; ++i_elem)
877 if (domain->elementsList->getElementType(i_elem) ==
static_cast<int>(tbox::ElType::HEX8))
879 max_element_size = 8;
882 if (domain->elementsList->getElementType(i_elem) ==
static_cast<int>(tbox::ElType::TETRA4))
884 max_element_size = 4;
889 K = Kokkos::View<scalar ***, Kokkos::LayoutRight, Kokkos::DefaultExecutionSpace>(
"A",
891 3 * max_element_size,
892 3 * max_element_size);
893 L = Kokkos::View<scalar ***, Kokkos::LayoutRight, Kokkos::DefaultExecutionSpace>(
"A",
897 S = Kokkos::View<scalar ***, Kokkos::LayoutRight, Kokkos::DefaultExecutionSpace>(
"A",
899 3 * max_element_size,
907 Kokkos::parallel_for(
908 numMyElems, KOKKOS_LAMBDA(
const int i_elem) {
909 bool bContinue =
true;
911 #ifdef KOKKOS_ENABLE_DEPRECATED_CODE
912 int i_thread = Kokkos::DefaultExecutionSpace::thread_pool_rank();
914 int i_thread = Kokkos::DefaultExecutionSpace::impl_thread_pool_rank();
916 const int e_size = domain->elementsList->getElementSize(i_elem);
918 auto K_e = subview(K, i_thread, make_pair(0, 3 * e_size), make_pair(0, 3 * e_size));
919 auto L_e = subview(L, i_thread, make_pair(0, e_size), make_pair(0, e_size));
920 auto S_e = subview(S, i_thread, make_pair(0, 3 * e_size), make_pair(0, e_size));
922 if (domain->elementsList->getElementType(i_elem) ==
static_cast<int>(tbox::ElType::HEX8))
924 elementMatrices_HEX8->compute(i_elem, i_thread);
925 Kokkos::deep_copy(K_e, elementMatrices_HEX8->getK(i_thread));
926 Kokkos::deep_copy(L_e, elementMatrices_HEX8->getL(i_thread));
927 Kokkos::deep_copy(S_e, elementMatrices_HEX8->getS(i_thread));
930 else if (domain->elementsList->getElementType(i_elem) ==
static_cast<int>(tbox::ElType::TETRA4))
932 elementMatrices_TETRA4->compute(i_elem, i_thread);
933 Kokkos::deep_copy(K_e, elementMatrices_TETRA4->getK(i_thread));
934 Kokkos::deep_copy(L_e, elementMatrices_TETRA4->getL(i_thread));
935 Kokkos::deep_copy(S_e, elementMatrices_TETRA4->getS(i_thread));
941 for (
auto i = 0; i < e_size; ++i)
944 global_ordinal_type global_node_id_i = algebraic->map->mapNodesWO->getGlobalElement(node_i);
946 if (algebraic->map->mapNodes->isNodeGlobalElement(global_node_id_i))
948 local_ordinal_type local_node_id_i = algebraic->map->mapNodes->getLocalElement(global_node_id_i);
949 for (
auto j = 0; j < e_size; ++j)
953 if (numPrimalDPN == 4 || numPrimalDPN == 1)
956 scalar vals = L_e(i, j);
958 local_L.sumIntoValues(local_node_id_i, &cols, 1, &vals,
false,
true);
960 if (numPrimalDPN == 4 || numPrimalDPN == 3)
962 for (
auto k = 0; k < 3; ++k)
963 for (
auto l = 0; l < 3; ++l)
966 scalar vals = K_e(3 * i + k, 3 * j + l);
967 local_K.sumIntoValues(3 * local_node_id_i + k, &cols, 1, &vals,
false,
true);
970 if (numPrimalDPN == 4)
972 for (
auto k = 0; k < 3; ++k)
975 scalar vals = S_e(3 * i + k, j);
976 local_S.sumIntoValues(3 * local_node_id_i + k, &cols, 1, &vals,
false,
true);
984 if (numPrimalDPN == 4 || numPrimalDPN == 1)
986 algebraic->matrices->
L->fillComplete();
988 if (numPrimalDPN == 4 || numPrimalDPN == 3)
990 algebraic->matrices->K->fillComplete();
992 if (numPrimalDPN == 4)
994 algebraic->matrices->S->fillComplete();
997 delete elementMatrices_HEX8;
998 delete elementMatrices_TETRA4;
1000 if (numPrimalDPN == 4)
1002 RCP<xpetra_bcrs_type> bA = rcp(
new xpetra_bcrs_type(algebraic->map->thermomecBlockedMap, algebraic->map->thermomecBlockedMap, 8));
1004 RCP<xpetra_crs_type> xL = rcp(
new xpetra_tcrs_type(algebraic->matrices->L));
1005 RCP<xpetra_crs_type> xK = rcp(
new xpetra_tcrs_type(algebraic->matrices->K));
1006 RCP<xpetra_crs_type> xS = rcp(
new xpetra_tcrs_type(algebraic->matrices->S));
1012 bA->setMatrix(0, 0, xwL);
1013 bA->setMatrix(1, 0, xwS);
1014 bA->setMatrix(1, 1, xwK);
1018 algebraic->matrices->A = bA;
1020 else if (numPrimalDPN == 3)
1022 RCP<xpetra_crs_type> xK = rcp(
new xpetra_tcrs_type(algebraic->matrices->K));
1027 RCP<xpetra_crs_type> xL = rcp(
new xpetra_tcrs_type(algebraic->matrices->L));
1040 template <
typename scalar>
1046 using Teuchos::tuple;
1048 typedef Kokkos::View<scalar ***, Kokkos::LayoutRight, Kokkos::DefaultExecutionSpace> elem_matrix_view_type;
1050 elem_matrix_view_type N, N_S;
1052 auto numMyElems = domain->elementsList->getElementNumber();
1054 #ifdef KOKKOS_ENABLE_DEPRECATED_CODE
1055 int pool_size = Kokkos::DefaultExecutionSpace::thread_pool_size();
1057 int pool_size = Kokkos::DefaultExecutionSpace::impl_thread_pool_size();
1059 int max_element_size = 1;
1060 int max_element_size_S = 1;
1062 for (
int i_elem = 0; i_elem < numMyElems; ++i_elem)
1064 if (domain->elementsList->getElementType(i_elem) ==
static_cast<int>(tbox::ElType::HEX8))
1066 max_element_size = 8;
1067 max_element_size_S = 4;
1070 if (domain->elementsList->getElementType(i_elem) ==
static_cast<int>(tbox::ElType::TETRA4))
1072 max_element_size = 4;
1073 max_element_size_S = 3;
1077 N = elem_matrix_view_type(
"A", pool_size, max_element_size, max_element_size);
1078 N_S = elem_matrix_view_type(
"A", pool_size, max_element_size_S, max_element_size_S);
1089 RCP<tpetra_vector_type> bWO = rcp(
new tpetra_vector_type(algebraic->map->mapDofsWO,
true));
1091 auto bWO_view = bWO->getLocalViewHost();
1099 auto numMySources = loads->sourcesList->getSourceNumber();
1101 for (
auto i = 0; i < numMySources; ++i)
1103 auto currentSourceSize = loads->sourcesList->getSourceSize(i);
1104 Kokkos::parallel_for(
1105 currentSourceSize, KOKKOS_LAMBDA(
const int e1) {
1107 #ifdef KOKKOS_ENABLE_DEPRECATED_CODE
1108 int i_thread = Kokkos::DefaultExecutionSpace::thread_pool_rank();
1110 int i_thread = Kokkos::DefaultExecutionSpace::impl_thread_pool_rank();
1112 const int e_size = domain->elementsList->getElementSize(i_elem);
1114 auto N_e = subview(N, i_thread, make_pair(0, e_size), make_pair(0, e_size));
1116 bool bContinue =
true;
1118 if (domain->elementsList->getElementType(i_elem) ==
static_cast<int>(tbox::ElType::HEX8))
1120 elementVectors_HEX8->compute(i_elem, i_thread);
1121 Kokkos::deep_copy(N_e, elementVectors_HEX8->getN(i_thread));
1123 else if (domain->elementsList->getElementType(i_elem) ==
static_cast<int>(tbox::ElType::TETRA4))
1125 elementVectors_TETRA4->compute(i_elem, i_thread);
1126 Kokkos::deep_copy(N_e, elementVectors_TETRA4->getN(i_thread));
1133 for (
auto ii = 0; ii < e_size; ++ii)
1139 if (algebraic->map->mapNodes->isNodeGlobalElement(global_i))
1143 for (
auto jj = 0; jj < e_size; ++jj)
1146 tmp *= loads->sourcesList->getSourceValue(i);
1148 Kokkos::atomic_fetch_add(&bWO_view(local_i, 0), tmp);
1161 auto numMynBCs = loads->neumannList->getNeumannNumber();
1163 const size_t numNodes = algebraic->map->mapNodes->getGlobalNumElements();
1165 const size_t numMyNodesWO = algebraic->map->mapNodesWO->getNodeNumElements();
1167 for (
auto i = 0; i < numMynBCs; ++i)
1169 auto currentBCSize = loads->neumannList->getNeumannSize(i);
1170 Kokkos::parallel_for(
1171 currentBCSize, KOKKOS_LAMBDA(
const int e1) {
1174 #ifdef KOKKOS_ENABLE_DEPRECATED_CODE
1175 int i_thread = Kokkos::DefaultExecutionSpace::thread_pool_rank();
1177 int i_thread = Kokkos::DefaultExecutionSpace::impl_thread_pool_rank();
1179 const int e_size = domain->elementsList->getElementSize(i_elem);
1181 auto N_e = subview(N_S, i_thread, make_pair(0, e_size), make_pair(0, e_size));
1183 bool bContinue =
true;
1185 if (domain->elementsList->getElementType(i_elem) ==
static_cast<int>(tbox::ElType::TRI3))
1187 elementVectors_TRI3->compute(i_elem, i_thread);
1188 Kokkos::deep_copy(N_e, elementVectors_TRI3->getN(i_thread));
1190 else if (domain->elementsList->getElementType(i_elem) ==
static_cast<int>(tbox::ElType::QUAD4))
1192 elementVectors_QUAD4->compute(i_elem, i_thread);
1193 Kokkos::deep_copy(N_e, elementVectors_QUAD4->getN(i_thread));
1200 for (
auto j = 0; j < numPrimalDPN; ++j)
1202 scalar which_dof = loads->neumannList->getNeumannDof(i, j);
1203 scalar value = loads->neumannList->getNeumannValue(i, j);
1204 if (which_dof == (scalar)1)
1206 for (
auto ii = 0; ii < e_size; ++ii)
1210 local_i = domain->elementsList->getElementNode(i_elem, ii);
1213 if (numPrimalDPN == 4)
1214 local_i = numMyNodesWO + 3 * domain->elementsList->getElementNode(i_elem, ii) + j;
1216 local_i = 3 * domain->elementsList->getElementNode(i_elem, ii) + j;
1221 for (
auto jj = 0; jj < e_size; ++jj)
1226 Kokkos::atomic_fetch_add(&bWO_view(local_i, 0), tmp);
1240 Tpetra::Export<> exportBWO(algebraic->map->mapDofsWO, algebraic->map->mapDofs);
1241 algebraic->vectors->b->doExport(*bWO, exportBWO, Tpetra::ADD);
1243 if (numPrimalDPN == 4 || numPrimalDPN == 3)
1244 algebraic->matrices->K->resumeFill();
1245 if (numPrimalDPN == 4 || numPrimalDPN == 1)
1246 algebraic->matrices->L->resumeFill();
1247 if (numPrimalDPN == 4)
1248 algebraic->matrices->S->resumeFill();
1250 auto numMyBCs = loads->dirichletList->getDirichletNumber();
1251 for (
auto i = 0; i < numMyBCs; ++i)
1253 auto currentBCSize = loads->dirichletList->getDirichletSize(i);
1254 for (
int i_node = 0; i_node < currentBCSize; ++i_node)
1258 if (algebraic->map->mapNodes->isNodeGlobalElement(n1))
1259 for (
auto j = 0; j < numPrimalDPN; ++j)
1261 scalar which_dof = loads->dirichletList->getDirichletDof(i, j);
1263 if (which_dof == (scalar)1)
1267 bool thermo =
false;
1272 global_index_r0 = n1;
1278 if (numPrimalDPN == 4)
1279 global_index_0 = numNodes;
1282 global_index_r0 = 3 * n1 + j;
1285 global_index = global_index_0 + global_index_r0;
1287 Teuchos::Array<global_ordinal_type> global_indices;
1291 size_t numColInds = algebraic->graph->L->getNumEntriesInGlobalRow(global_index);
1292 global_indices.resize(numColInds);
1293 algebraic->graph->L->getGlobalRowCopy(global_index, global_indices(), numColInds);
1296 for (
auto k = 0; k < numColInds; ++k)
1297 algebraic->matrices->L->replaceGlobalValues(global_index,
1298 tuple<global_ordinal_type>(global_indices[k]),
1299 tuple<scalar>(scalar(0.0)));
1302 algebraic->matrices->L->replaceGlobalValues(global_index,
1303 tuple<global_ordinal_type>(global_index),
1304 tuple<scalar>(scalar(1.0)));
1308 size_t numColInds = algebraic->graph->K->getNumEntriesInGlobalRow(global_index_r0);
1309 global_indices.resize(numColInds);
1310 algebraic->graph->K->getGlobalRowCopy(global_index_r0, global_indices(), numColInds);
1313 for (
auto k = 0; k < numColInds; ++k)
1314 algebraic->matrices->K->replaceGlobalValues(global_index_r0,
1315 tuple<global_ordinal_type>(global_indices[k]),
1316 tuple<scalar>(scalar(0.0)));
1319 algebraic->matrices->K->replaceGlobalValues(global_index_r0,
1320 tuple<global_ordinal_type>(global_index_r0),
1321 tuple<scalar>(scalar(1.0)));
1323 if (numPrimalDPN == 4)
1325 size_t numColInds = algebraic->graph->S->getNumEntriesInGlobalRow(global_index_r0);
1326 global_indices.resize(numColInds);
1327 algebraic->graph->S->getGlobalRowCopy(global_index_r0, global_indices(), numColInds);
1330 for (
auto k = 0; k < numColInds; ++k)
1331 algebraic->matrices->S->replaceGlobalValues(global_index_r0,
1332 tuple<global_ordinal_type>(global_indices[k]),
1333 tuple<scalar>(scalar(0.0)));
1338 algebraic->vectors->b->replaceGlobalValue(global_index,
1339 (scalar(loads->dirichletList->getDirichletValue(i, j))));
1340 algebraic->vectors->x->replaceGlobalValue(global_index,
1341 (scalar(loads->dirichletList->getDirichletValue(i, j))));
1346 if (numPrimalDPN == 4 || numPrimalDPN == 3)
1347 algebraic->matrices->K->fillComplete();
1348 if (numPrimalDPN == 4 || numPrimalDPN == 1)
1349 algebraic->matrices->L->fillComplete();
1350 if (numPrimalDPN == 4)
1351 algebraic->matrices->S->fillComplete();
1381 auto numMyWeightRegions = loads->weightsList->getWeightRegionsNumber();
1383 for (
auto i = 0; i < numMyWeightRegions; ++i)
1385 auto currentWeightRegionSize = loads->weightsList->getWeightRegionSize(i);
1386 Kokkos::parallel_for(
1387 currentWeightRegionSize, KOKKOS_LAMBDA(
const int i_node) {
1390 if (algebraic->map->mapNodes->isNodeGlobalElement(n1))
1392 for (auto j = 0; j < numPrimalDPN; ++j)
1394 if (loads->weightsList->getWeightDof(i, j))
1396 global_ordinal_type global_index;
1404 if (numPrimalDPN == 4)
1405 global_index = numNodes + 3 * n1 + j;
1407 global_index = 3 * n1 + j;
1410 algebraic->vectors->weights->sumIntoGlobalValue(global_index,
1411 (scalar(loads->weightsList->getWeightValue(i, j))),
1419 delete elementVectors_TRI3;
1420 delete elementVectors_QUAD4;
1421 delete elementVectors_TETRA4;
1422 delete elementVectors_HEX8;