waves
Basic FE playground
DiscreteProblem.hpp
Go to the documentation of this file.
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>
8 
9 namespace katoptron
10 {
11 
27 template <typename scalar>
29  Teuchos::RCP<Teuchos::ParameterList> randomParams,
30  Kokkos::View<scalar *, Kokkos::LayoutLeft> m_rv_values)
31 {
32  numPrimalDPN = _numPrimalDPN;
33 
34  LinearSolver::getTimers()["map"].start();
35  Teuchos::RCP<Map> map = Teuchos::rcp(new Map(pbl, numPrimalDPN));
36  LinearSolver::getTimers()["map"].stop();
37  LinearSolver::getTimers()["domain"].start();
38  domain = Teuchos::rcp(new Domain<scalar>(pbl, map, randomParams, m_rv_values));
39  LinearSolver::getTimers()["domain"].stop();
40  LinearSolver::getTimers()["loads"].start();
41  loads = Teuchos::rcp(new Loads<scalar>(pbl, map, domain));
42  LinearSolver::getTimers()["loads"].stop();
43  algebraic = Teuchos::rcp(new Algebraic<scalar>(map, domain));
44 }
45 
55 template <typename scalar>
56 Mask<scalar> DiscreteProblem<scalar>::updateActivity(int active_set_iteration)
57 {
58  using Teuchos::Array;
59  using Teuchos::RCP;
60  using Teuchos::rcp;
61 
62  const size_t myRank = algebraic->map->mapDofs->getComm()->getRank();
63 
64  // First, get the displacement solution,
65 
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);
69 
70  RCP<tpetra_vector_type> x = rcp(new tpetra_vector_type(algebraic->map->mapDofs, false));
71 
72  RCP<xpetra_bmvector_type> blockedSol =
73  rcp(new xpetra_bmvector_type(algebraic->map->blockedMap, algebraic->vectors->solutionMultiVector));
74 
75  RCP<tpetra_vector_type> l =
76  (Teuchos::rcp_dynamic_cast<xpetra_tmvector_type>(blockedSol->getMultiVector(1, true)))->getTpetra_MultiVector()->getVectorNonConst(0);
77 
78  RCP<tpetra_vector_type> gap = rcp(new tpetra_vector_type(*algebraic->vectors->initialGap, Teuchos::Copy));
79 
80  RCP<export_type> exportx = rcp(new export_type(algebraic->map->fullmap, algebraic->map->mapDofs));
81 
82  x->doExport(*xl, *exportx, Tpetra::INSERT);
83 
84  // Multiply it with B_G and remove initial gap
85 
86  // gap = 1 * B_G * x - 1 * gap
87 
88  algebraic->matrices->B_G->apply(*x, *gap, Teuchos::NO_TRANS, ((scalar)1.), ((scalar)-1.));
89 
90  // sum lambda + c * the previously computed vector
91  double c = 1.;
92  l->update((scalar)c, *gap, (scalar)1.); // l = l + c*gap
93 
94  // loop on the local lagrange and create the new activity
95 
96  l->template sync<Kokkos::HostSpace>();
97  auto l_2d = l->template getLocalView<Kokkos::HostSpace>();
98 
99  Mask<scalar> hasConverged;
100  const size_t ensemble_size = hasConverged.getSize();
101 
102  for (size_t l = 0; l < ensemble_size; ++l)
103  hasConverged.set(l, true);
104 
105  bool hasConverged_old = false;
106  if (active_set_iteration > 1)
107  hasConverged_old = true;
108 
109  algebraic->matrices->B_2->resumeFill();
110  algebraic->matrices->C->resumeFill();
111  algebraic->matrices->Cb->resumeFill();
112 
113  const size_t numDPN = algebraic->map->numPrimalDPN;
114 
115  size_t numLMPN;
116  if (domain->contactsList->hasAtLeastOneSticking() && numDPN == 4)
117  numLMPN = 4;
118  else if (domain->contactsList->hasAtLeastOneSticking())
119  numLMPN = 3;
120  else if (numDPN == 4)
121  numLMPN = 2;
122  else
123  numLMPN = 1;
124 
125  const int numMyLagrangeNum = algebraic->map->mapLagrangeDofs->getNodeNumElements();
126 
127  const int numMyMechLagrangeNum = numMyLagrangeNum / numLMPN;
128 
129  for (auto local_index = 0; local_index < numMyMechLagrangeNum; ++local_index)
130  {
131  global_ordinal_type lagrange_id = algebraic->map->mapLagrangeDofs->getGlobalElement(local_index);
132  global_ordinal_type node_id = algebraic->map->lm_to_dof[lagrange_id];
133  lagrange_id = lagrange_id * numLMPN;
134  local_ordinal_type interface_id = domain->contactsList->getInterfaceOfSlaveNode(node_id);
135 
136  if (domain->contactsList->isNotUpdated(interface_id))
137  continue;
138 
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.);
142 
143  scalar oneMask, oneNotMask;
144  mask_assign<scalar>(current_mask, oneMask) = {1., 0.};
145  mask_assign<scalar>(current_mask, oneNotMask) = {0., 1.};
146 
147  if (active_set_iteration > 1)
148  if (!MaskLogic::AND(oneMask == old_old_mask))
149  hasConverged_old = false;
150 
151  hasConverged = hasConverged && (oneMask == old_mask);
152 
153  if (!MaskLogic::AND(oneMask == old_mask))
154  {
155  global_ordinal_type lagrange_id = algebraic->map->mapLagrangeDofs->getGlobalElement(local_index);
156 
157  lagrange_id = lagrange_id * numLMPN;
158 
159  // If the current Lagrange multiplier is active for every sample of the current ensemble,
160  // we can do the update without using the mask entries.
161  if (MaskLogic::AND(current_mask))
162  {
163  Teuchos::ArrayView<const local_ordinal_type> local_indices;
164  Teuchos::ArrayView<const scalar> values;
165 
166  if (domain->contactsList->isSticking(interface_id))
167  {
168  for (size_t j = 0; j < numLMPN; ++j)
169  {
170  algebraic->matrices->B->getLocalRowView(local_index * numLMPN + j, local_indices, values);
171 
172  for (auto k = 0; k < local_indices.size(); ++k)
173  {
174  global_ordinal_type column_index = algebraic->matrices->B->getColMap()->getGlobalElement(local_indices[k]);
175 
176  algebraic->matrices->B_2->replaceGlobalValues(lagrange_id + j,
177  Teuchos::tuple<global_ordinal_type>(column_index),
178  Teuchos::tuple<scalar>(values[k]));
179  }
180 
181  algebraic->matrices->C->replaceGlobalValues(lagrange_id + j,
182  Teuchos::tuple<global_ordinal_type>(lagrange_id + j),
183  Teuchos::tuple<scalar>((scalar)0.));
184 
185  algebraic->matrices->Cb->replaceGlobalValues(lagrange_id + j,
186  Teuchos::tuple<global_ordinal_type>(lagrange_id + j),
187  Teuchos::tuple<scalar>((scalar)1.));
188  }
189  }
190  else
191  {
192  algebraic->matrices->B->getLocalRowView(local_index * numLMPN, local_indices, values);
193 
194  for (auto k = 0; k < local_indices.size(); ++k)
195  {
196  global_ordinal_type column_index = algebraic->matrices->B->getColMap()->getGlobalElement(local_indices[k]);
197 
198  algebraic->matrices->B_2->replaceGlobalValues(lagrange_id,
199  Teuchos::tuple<global_ordinal_type>(column_index),
200  Teuchos::tuple<scalar>(values[k]));
201  }
202 
203  algebraic->matrices->C->replaceGlobalValues(lagrange_id,
204  Teuchos::tuple<global_ordinal_type>(lagrange_id),
205  Teuchos::tuple<scalar>((scalar)0.));
206 
207  algebraic->matrices->Cb->replaceGlobalValues(lagrange_id,
208  Teuchos::tuple<global_ordinal_type>(lagrange_id),
209  Teuchos::tuple<scalar>((scalar)1.));
210 
211  if (numLMPN == 2 || numLMPN == 4)
212  {
213  algebraic->matrices->B->getLocalRowView(local_index * numLMPN + numLMPN - 1, local_indices, values);
214 
215  for (auto k = 0; k < local_indices.size(); ++k)
216  {
217  global_ordinal_type column_index = algebraic->matrices->B->getColMap()->getGlobalElement(local_indices[k]);
218 
219  algebraic->matrices->B_2->replaceGlobalValues(lagrange_id + numLMPN - 1,
220  Teuchos::tuple<global_ordinal_type>(column_index),
221  Teuchos::tuple<scalar>(values[k]));
222  }
223 
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.));
227 
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.));
231  }
232  }
233  }
234  // If the current Lagrange multiplier is inactive for every sample of the current ensemble,
235  // we can do the update without using the mask entries.
236  else if (!MaskLogic::AND(current_mask))
237  {
238  Teuchos::ArrayView<const local_ordinal_type> local_indices;
239  Teuchos::ArrayView<const scalar> values;
240 
241  if (domain->contactsList->isSticking(interface_id))
242  {
243  for (size_t j = 0; j < numLMPN; ++j)
244  {
245  algebraic->matrices->B->getLocalRowView(local_index * numLMPN + j, local_indices, values);
246 
247  for (auto k = 0; k < local_indices.size(); ++k)
248  {
249  global_ordinal_type column_index = algebraic->matrices->B->getColMap()->getGlobalElement(local_indices[k]);
250 
251  algebraic->matrices->B_2->replaceGlobalValues(lagrange_id + j,
252  Teuchos::tuple<global_ordinal_type>(column_index),
253  Teuchos::tuple<scalar>((scalar)0.));
254  }
255 
256  algebraic->matrices->C->replaceGlobalValues(lagrange_id + j,
257  Teuchos::tuple<global_ordinal_type>(lagrange_id + j),
258  Teuchos::tuple<scalar>((scalar)1.));
259 
260  algebraic->matrices->Cb->replaceGlobalValues(lagrange_id + j,
261  Teuchos::tuple<global_ordinal_type>(lagrange_id + j),
262  Teuchos::tuple<scalar>((scalar)0.));
263  }
264  }
265  else
266  {
267  algebraic->matrices->B->getLocalRowView(local_index * numLMPN, local_indices, values);
268 
269  for (auto k = 0; k < local_indices.size(); ++k)
270  {
271  global_ordinal_type column_index = algebraic->matrices->B->getColMap()->getGlobalElement(local_indices[k]);
272 
273  algebraic->matrices->B_2->replaceGlobalValues(lagrange_id,
274  Teuchos::tuple<global_ordinal_type>(column_index),
275  Teuchos::tuple<scalar>((scalar)0.));
276  }
277 
278  algebraic->matrices->C->replaceGlobalValues(lagrange_id,
279  Teuchos::tuple<global_ordinal_type>(lagrange_id),
280  Teuchos::tuple<scalar>((scalar)1.));
281 
282  algebraic->matrices->Cb->replaceGlobalValues(lagrange_id,
283  Teuchos::tuple<global_ordinal_type>(lagrange_id),
284  Teuchos::tuple<scalar>((scalar)0.));
285 
286  if (numLMPN == 2 || numLMPN == 4)
287  {
288  algebraic->matrices->B->getLocalRowView(local_index * numLMPN + numLMPN - 1, local_indices, values);
289 
290  for (auto k = 0; k < local_indices.size(); ++k)
291  {
292  global_ordinal_type column_index = algebraic->matrices->B->getColMap()->getGlobalElement(local_indices[k]);
293 
294  algebraic->matrices->B_2->replaceGlobalValues(lagrange_id + numLMPN - 1,
295  Teuchos::tuple<global_ordinal_type>(column_index),
296  Teuchos::tuple<scalar>((scalar)0.));
297  }
298 
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.));
302 
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.));
306  }
307  }
308  }
309  // Else, if the activity of the current Lagrange multiplier depends on the samples, we use the mask to do the update.
310  else
311  {
312 
313  Teuchos::ArrayView<const local_ordinal_type> local_indices;
314  Teuchos::ArrayView<const scalar> values;
315 
316  if (domain->contactsList->isSticking(interface_id))
317  {
318  for (size_t j = 0; j < numLMPN; ++j)
319  {
320  algebraic->matrices->B->getLocalRowView(local_index * numLMPN + j, local_indices, values);
321 
322  for (auto k = 0; k < local_indices.size(); ++k)
323  {
324  global_ordinal_type column_index = algebraic->matrices->B->getColMap()->getGlobalElement(local_indices[k]);
325 
326  algebraic->matrices->B_2->replaceGlobalValues(lagrange_id + j,
327  Teuchos::tuple<global_ordinal_type>(column_index),
328  Teuchos::tuple<scalar>(oneMask * values[k]));
329  }
330 
331  algebraic->matrices->C->replaceGlobalValues(lagrange_id + j,
332  Teuchos::tuple<global_ordinal_type>(lagrange_id + j),
333  Teuchos::tuple<scalar>(oneNotMask));
334 
335  algebraic->matrices->Cb->replaceGlobalValues(lagrange_id + j,
336  Teuchos::tuple<global_ordinal_type>(lagrange_id + j),
337  Teuchos::tuple<scalar>(oneMask));
338  }
339  }
340  else
341  {
342  algebraic->matrices->B->getLocalRowView(local_index * numLMPN, local_indices, values);
343 
344  for (auto k = 0; k < local_indices.size(); ++k)
345  {
346  global_ordinal_type column_index = algebraic->matrices->B->getColMap()->getGlobalElement(local_indices[k]);
347 
348  algebraic->matrices->B_2->replaceGlobalValues(lagrange_id,
349  Teuchos::tuple<global_ordinal_type>(column_index),
350  Teuchos::tuple<scalar>(oneMask * values[k]));
351  }
352 
353  algebraic->matrices->C->replaceGlobalValues(lagrange_id,
354  Teuchos::tuple<global_ordinal_type>(lagrange_id),
355  Teuchos::tuple<scalar>(oneNotMask));
356 
357  algebraic->matrices->Cb->replaceGlobalValues(lagrange_id,
358  Teuchos::tuple<global_ordinal_type>(lagrange_id),
359  Teuchos::tuple<scalar>(oneMask));
360 
361  if (numLMPN == 2 || numLMPN == 4)
362  {
363  algebraic->matrices->B->getLocalRowView(local_index * numLMPN + numLMPN - 1, local_indices, values);
364 
365  for (auto k = 0; k < local_indices.size(); ++k)
366  {
367  global_ordinal_type column_index = algebraic->matrices->B->getColMap()->getGlobalElement(local_indices[k]);
368 
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]));
372  }
373 
374  algebraic->matrices->C->replaceGlobalValues(lagrange_id + numLMPN - 1,
375  Teuchos::tuple<global_ordinal_type>(lagrange_id + numLMPN - 1),
376  Teuchos::tuple<scalar>(oneNotMask));
377 
378  algebraic->matrices->Cb->replaceGlobalValues(lagrange_id + numLMPN - 1,
379  Teuchos::tuple<global_ordinal_type>(lagrange_id + numLMPN - 1),
380  Teuchos::tuple<scalar>(oneMask));
381  }
382  }
383  }
384  }
385 
386  old_old_activity(local_index) = old_activity(local_index);
387  old_activity(local_index) = oneMask;
388  }
389 
390  algebraic->matrices->B_2->fillComplete(algebraic->map->mapDofs, algebraic->map->mapLagrangeDofs);
391  algebraic->matrices->C->fillComplete();
392  algebraic->matrices->Cb->fillComplete();
393 
394  RCP<tpetra_vector_type> initialActiveGap(new tpetra_vector_type(algebraic->map->mapLagrangeDofs, true));
395  algebraic->matrices->Cb->apply(*algebraic->vectors->initialGap, *initialActiveGap);
396 
397  RCP<xpetra_mvector_type> xg =
398  rcp(new xpetra_tmvector_type(initialActiveGap));
399  algebraic->vectors->rhsBlockedMultiVector->setMultiVector(1, xg, true);
400  algebraic->vectors->rhsMultiVector = algebraic->vectors->rhsBlockedMultiVector->Merge();
401 
402  size_t numtasks = algebraic->map->mapDofs->getComm()->getSize();
403 
404  double tmp;
405 
406  if (myRank == 0)
407  {
408  Mask<scalar> globalHasConverged = hasConverged;
409  bool globalHasConverged_old = hasConverged_old;
410 
411  for (size_t i = 1; i < numtasks; ++i)
412  {
413  for (size_t l = 0; l < ensemble_size; ++l)
414  {
415  MPI_Recv(&tmp, 1, MPI_DOUBLE, i, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
416  if (tmp < 0.)
417  globalHasConverged.set(l, false);
418  }
419  MPI_Recv(&tmp, 1, MPI_DOUBLE, i, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
420  if (tmp < 0.)
421  globalHasConverged_old = false;
422  }
423  hasConverged = globalHasConverged;
424  hasConverged_old = globalHasConverged_old;
425 
426  for (size_t i = 1; i < numtasks; ++i)
427  {
428  for (size_t l = 0; l < ensemble_size; ++l)
429  {
430  tmp = hasConverged.get(l) ? 1. : -1.;
431  MPI_Send(&tmp, 1, MPI_DOUBLE, i, 0, MPI_COMM_WORLD);
432  }
433  tmp = hasConverged_old ? 1. : -1.;
434  MPI_Send(&tmp, 1, MPI_DOUBLE, i, 0, MPI_COMM_WORLD);
435  }
436  }
437  else
438  {
439  for (size_t l = 0; l < ensemble_size; ++l)
440  {
441  tmp = hasConverged.get(l) ? 1. : -1.;
442  MPI_Send(&tmp, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
443  }
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)
447  {
448  MPI_Recv(&tmp, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
449  hasConverged.set(l, tmp > 0.);
450  }
451  MPI_Recv(&tmp, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
452  hasConverged_old = (tmp > 0.);
453  }
454 
455  if (hasConverged_old)
456  {
457  if (myRank == 0)
458  {
459  std::cout << "-----------------------" << std::endl;
460  std::cout << "---Has not converged---" << std::endl;
461  std::cout << "-----------------------" << std::endl;
462  }
463  return true;
464  }
465 
466  return hasConverged;
467 }
468 
478 template <typename scalar>
479 void DiscreteProblem<scalar>::createBlockMatrix(scalar gamma, bool write, bool merge)
480 {
481  using Teuchos::Array;
482  using Teuchos::RCP;
483  using Teuchos::rcp;
484 
485  const std::vector<RCP<const tpetra_map_type>> maps = {algebraic->map->mapDofs, algebraic->map->mapLagrangeDofs};
486 
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();
492 
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();
496 
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];
502 
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];
508 
509  algebraic->map->fullmap =
510  rcp(new tpetra_map_type(numDofs + numLagrangeDofs,
511  myFullDofs,
512  algebraic->map->mapDofs->getIndexBase(),
513  algebraic->map->mapDofs->getComm()));
514  algebraic->map->fullmapWO =
515  rcp(new tpetra_map_type(numDofs + numLagrangeDofs,
516  myFullDofsWO,
517  algebraic->map->mapDofs->getIndexBase(),
518  algebraic->map->mapDofs->getComm()));
519 
520  RCP<xpetra_map_type> xfullmap = rcp(new xpetra_tmap_type(algebraic->map->fullmap));
521 
522  const size_t myRank = algebraic->map->mapDofs->getComm()->getRank();
523  const size_t numLclIndices = (myRank == 0) ? numDofs + numLagrangeDofs : 0;
524 
525  algebraic->map->fullmapOutput =
526  rcp(new tpetra_map_type(numDofs + numLagrangeDofs,
527  numLclIndices,
528  algebraic->map->mapDofs->getIndexBase(),
529  algebraic->map->mapDofs->getComm()));
530 
531  std::vector<RCP<const xpetra_map_type>> xsubmaps = {rcp(new xpetra_tmap_type(algebraic->map->mapDofs)), //algebraic->map->thermomecBlockedMap,
532  rcp(new xpetra_tmap_type(algebraic->map->mapLagrangeDofs))};
533 
534  algebraic->map->blockedMap =
535  rcp(new xpetra_bmap_type(xfullmap, xsubmaps, true));
536 
537  algebraic->matrices->blockedMatrix =
538  rcp(new xpetra_bcrs_type(algebraic->map->blockedMap, algebraic->map->blockedMap, 8));
539 
540  //Update B based on the initial activity
541 
542  algebraic->matrices->B_2->resumeFill();
543  algebraic->matrices->C->resumeFill();
544  algebraic->matrices->Cb->resumeFill();
545 
546  const size_t numDPN = algebraic->map->numPrimalDPN;
547 
548  size_t numLMPN;
549  if (domain->contactsList->hasAtLeastOneSticking() && numDPN == 4)
550  numLMPN = 4;
551  else if (domain->contactsList->hasAtLeastOneSticking())
552  numLMPN = 3;
553  else if (numDPN == 4)
554  numLMPN = 2;
555  else
556  numLMPN = 1;
557 
558  const int numMyLagrangeNum = algebraic->map->mapLagrangeDofs->getNodeNumElements();
559 
560  const int numMyMechLagrangeNum = numMyLagrangeNum / numLMPN;
561 
562  old_activity = Kokkos::View<scalar *, Kokkos::LayoutRight>("R", numMyMechLagrangeNum);
563  old_old_activity = Kokkos::View<scalar *, Kokkos::LayoutRight>("R", numMyMechLagrangeNum);
564 
565  for (auto local_index = 0; local_index < numMyMechLagrangeNum; ++local_index)
566  {
567  global_ordinal_type lagrange_id = algebraic->map->mapLagrangeDofs->getGlobalElement(local_index);
568  global_ordinal_type node_id = algebraic->map->lm_to_dof[lagrange_id];
569  lagrange_id = lagrange_id * numLMPN;
570  local_ordinal_type interface_id = domain->contactsList->getInterfaceOfSlaveNode(node_id);
571  if (domain->contactsList->isInitiallyOpen(interface_id))
572  {
573  if (domain->contactsList->isNodeInitiallyClosed(interface_id, node_id))
574  {
575  old_activity(local_index) = (scalar)1.;
576  continue;
577  }
578  old_activity(local_index) = (scalar)0.;
579 
580  Teuchos::ArrayView<const local_ordinal_type> local_indices;
581  Teuchos::ArrayView<const scalar> values;
582 
583  for (size_t j = 0; j < numLMPN; ++j)
584  {
585  algebraic->matrices->B_2->getLocalRowView(local_index * numLMPN + j, local_indices, values);
586 
587  for (auto k = 0; k < local_indices.size(); ++k)
588  {
589  global_ordinal_type column_index = algebraic->matrices->B->getColMap()->getGlobalElement(local_indices[k]);
590 
591  algebraic->matrices->B_2->replaceGlobalValues(lagrange_id + j,
592  Teuchos::tuple<global_ordinal_type>(column_index),
593  Teuchos::tuple<scalar>((scalar)0.));
594  }
595 
596  algebraic->matrices->C->replaceGlobalValues(lagrange_id + j,
597  Teuchos::tuple<global_ordinal_type>(lagrange_id + j),
598  Teuchos::tuple<scalar>((scalar)1.));
599 
600  algebraic->matrices->Cb->replaceGlobalValues(lagrange_id + j,
601  Teuchos::tuple<global_ordinal_type>(lagrange_id + j),
602  Teuchos::tuple<scalar>((scalar)0.));
603  }
604  }
605  else
606  {
607  old_activity(local_index) = (scalar)1.;
608  if (domain->contactsList->isSticking(interface_id) == false && domain->contactsList->hasAtLeastOneSticking() == true)
609  {
610  Teuchos::ArrayView<const local_ordinal_type> local_indices;
611  Teuchos::ArrayView<const scalar> values;
612 
613  for (size_t j = 1; j < numLMPN; ++j)
614  {
615  algebraic->matrices->B_2->getLocalRowView(local_index * numLMPN + j, local_indices, values);
616 
617  for (auto k = 0; k < local_indices.size(); ++k)
618  {
619  global_ordinal_type column_index = algebraic->matrices->B->getColMap()->getGlobalElement(local_indices[k]);
620 
621  algebraic->matrices->B_2->replaceGlobalValues(lagrange_id + j,
622  Teuchos::tuple<global_ordinal_type>(column_index),
623  Teuchos::tuple<scalar>((scalar)0.));
624  }
625 
626  algebraic->matrices->C->replaceGlobalValues(lagrange_id + j,
627  Teuchos::tuple<global_ordinal_type>(lagrange_id + j),
628  Teuchos::tuple<scalar>((scalar)1.));
629 
630  algebraic->matrices->Cb->replaceGlobalValues(lagrange_id + j,
631  Teuchos::tuple<global_ordinal_type>(lagrange_id + j),
632  Teuchos::tuple<scalar>((scalar)0.));
633  }
634  }
635  }
636  }
637  algebraic->matrices->B_2->fillComplete(algebraic->map->mapDofs, algebraic->map->mapLagrangeDofs);
638  algebraic->matrices->C->fillComplete();
639  algebraic->matrices->Cb->fillComplete();
640 
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));
644 
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));
648 
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));
654  CTC =
655  Xpetra::MatrixMatrix<scalar, local_ordinal_type, global_ordinal_type, node_type>::Multiply(*(algebraic->matrices->xwB_T),
656  false,
657  *(algebraic->matrices->xwB),
658  false,
659  *verbOut);
660  RCP<xpetra_bcrs_type> bA = Teuchos::rcp_dynamic_cast<xpetra_bcrs_type>(algebraic->matrices->A);
661 
662  if (gamma > 0.)
663  {
664  if (bA != Teuchos::null && merge)
665  Xpetra::MatrixMatrix<scalar, local_ordinal_type, global_ordinal_type, node_type>::TwoMatrixAdd(*CTC,
666  false,
667  gamma,
668  *(bA->Merge()),
669  false,
670  (scalar)1.,
671  ApCTC,
672  *verbOut);
673  else
674  Xpetra::MatrixMatrix<scalar, local_ordinal_type, global_ordinal_type, node_type>::TwoMatrixAdd(*CTC,
675  false,
676  gamma,
677  *(algebraic->matrices->A),
678  false,
679  (scalar)1.,
680  ApCTC,
681  *verbOut);
682  }
683  else
684  {
685  if (bA != Teuchos::null && merge)
686  ApCTC = bA->Merge();
687  else
688  ApCTC = algebraic->matrices->A;
689  }
690 
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);
695 
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();
700 
701  if (write)
702  Tpetra::MatrixMarket::Writer<tpetra_crs_type>::writeSparseFile(std::string("ApCTC_mm.txt"), tpetra_ApCTC);
703 }
704 
713 template <typename scalar>
715 {
716  using Teuchos::Array;
717  using Teuchos::RCP;
718  using Teuchos::rcp;
719 
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));
723 
724  algebraic->vectors->lagrange = rcp(new tpetra_vector_type(algebraic->map->mapLagrangeDofs, true));
725 
726  {
727  size_t numMyContacts = domain->contactsList->getContactNumber();
728  size_t numPrimalDPN = algebraic->map->numPrimalDPN;
729  size_t numMPrimalDPN = 3;
730  size_t offSet = 0;
731  if (numPrimalDPN == 4)
732  offSet = algebraic->map->mapNodes->getGlobalNumElements();
733 
734  for (auto i = 0; i < numMyContacts; ++i)
735  {
736  size_t numMySlaveNodes = domain->contactsList->getSlaveNodesSize(i);
737  for (auto j = 0; j < numMySlaveNodes; ++j)
738  {
739  // get unique node id
740  auto global_nodeid = domain->contactsList->getSlaveNode(i, j);
741  if (algebraic->map->mapNodes->isNodeGlobalElement(global_nodeid))
742  {
743  auto local_nodeid = algebraic->map->mapNodesWO->getLocalElement(global_nodeid);
744  for (auto k = 0; k < 3; ++k)
745  {
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);
750  }
751  }
752  }
753 
754  size_t numMyMasterNodes = domain->contactsList->getMasterNodesSize(i);
755  for (auto j = 0; j < numMyMasterNodes; ++j)
756  {
757  // get unique node id
758  auto global_nodeid = domain->contactsList->getMasterNode(i, j);
759  if (algebraic->map->mapNodes->isNodeGlobalElement(global_nodeid))
760  {
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)));
765  }
766  }
767  }
768  }
769 
770  algebraic->matrices->B_G->apply(*intialContactPosition, *algebraic->vectors->initialGap, Teuchos::NO_TRANS, ((scalar)-1.));
771 
772  algebraic->matrices->Cb->apply(*algebraic->vectors->initialGap, *initialActiveGap);
773 
774  // TODO: fix me when there is mesh tying and contact together!
775  if (domain->contactsList->isTying(0))
776  initialActiveGap->scale((scalar)0.);
777 
778  if (write)
779  Tpetra::MatrixMarket::Writer<tpetra_crs_type>::writeDenseFile(std::string("initialActiveGap_mm.txt"), initialActiveGap);
780 
781  RCP<xpetra_mvector_type> xb =
782  rcp(new xpetra_tmvector_type(algebraic->vectors->b));
783 
784  RCP<xpetra_mvector_type> xg =
785  rcp(new xpetra_tmvector_type(initialActiveGap));
786 
787  RCP<xpetra_mvector_type> xx =
788  rcp(new xpetra_tmvector_type(algebraic->vectors->x));
789 
790  RCP<xpetra_mvector_type> xlagrange =
791  rcp(new xpetra_tmvector_type(algebraic->vectors->lagrange));
792 
793  algebraic->vectors->rhsBlockedMultiVector =
794  rcp(new xpetra_bmvector_type(algebraic->map->blockedMap, 1, true));
795 
796  algebraic->vectors->solutionBlockedMultiVector =
797  rcp(new xpetra_bmvector_type(algebraic->map->blockedMap, 1, true));
798 
799  algebraic->vectors->rhsBlockedMultiVector->setMultiVector(0, xb, true);
800  algebraic->vectors->rhsBlockedMultiVector->setMultiVector(1, xg, true);
801 
802  algebraic->vectors->solutionBlockedMultiVector->setMultiVector(0, xx, true);
803  algebraic->vectors->solutionBlockedMultiVector->setMultiVector(1, xlagrange, true);
804 
805  algebraic->vectors->rhsMultiVector = algebraic->vectors->rhsBlockedMultiVector->Merge();
806  algebraic->vectors->solutionMultiVector = algebraic->vectors->solutionBlockedMultiVector->Merge();
807 }
808 
817 template <typename scalar>
819 {
820  compute_B<scalar>(algebraic->matrices,
821  domain->elementsList,
822  domain->nodesList,
823  domain->contactsList,
824  algebraic->map, comm);
825 }
826 
830 template <typename scalar>
832 {
833  using Teuchos::RCP;
834  using Teuchos::rcp;
835 
836  typedef KokkosSparse::CrsMatrix<scalar, local_ordinal_type, typename tpetra_crs_type::execution_space> local_matrix_type;
837 
838  typedef Kokkos::View<scalar ***, Kokkos::LayoutRight, Kokkos::DefaultExecutionSpace> elem_matrix_view_type;
839 
840  local_matrix_type local_L, local_K, local_S;
841 
842  if (numPrimalDPN == 4 || numPrimalDPN == 1)
843  {
844  algebraic->matrices->L->resumeFill();
845  algebraic->matrices->L->setAllToScalar((scalar)0.0);
846 
847  local_L = algebraic->matrices->L->getLocalMatrix();
848  }
849  if (numPrimalDPN == 4 || numPrimalDPN == 3)
850  {
851  algebraic->matrices->K->resumeFill();
852  algebraic->matrices->K->setAllToScalar((scalar)0.0);
853 
854  local_K = algebraic->matrices->K->getLocalMatrix();
855  }
856  if (numPrimalDPN == 4)
857  {
858  algebraic->matrices->S->resumeFill();
859  algebraic->matrices->S->setAllToScalar((scalar)0.0);
860 
861  local_S = algebraic->matrices->S->getLocalMatrix();
862  }
863 
864  auto numMyElems = domain->elementsList->getElementNumber();
865 
866  elem_matrix_view_type K, L, S;
867 
868 #ifdef KOKKOS_ENABLE_DEPRECATED_CODE
869  int pool_size = Kokkos::DefaultExecutionSpace::thread_pool_size();
870 #else
871  int pool_size = Kokkos::DefaultExecutionSpace::impl_thread_pool_size();
872 #endif
873  int max_element_size = 1;
874 
875  for (int i_elem = 0; i_elem < numMyElems; ++i_elem)
876  {
877  if (domain->elementsList->getElementType(i_elem) == static_cast<int>(tbox::ElType::HEX8))
878  {
879  max_element_size = 8;
880  break;
881  }
882  if (domain->elementsList->getElementType(i_elem) == static_cast<int>(tbox::ElType::TETRA4))
883  {
884  max_element_size = 4;
885  break;
886  }
887  }
888 
889  K = Kokkos::View<scalar ***, Kokkos::LayoutRight, Kokkos::DefaultExecutionSpace>("A",
890  pool_size,
891  3 * max_element_size,
892  3 * max_element_size);
893  L = Kokkos::View<scalar ***, Kokkos::LayoutRight, Kokkos::DefaultExecutionSpace>("A",
894  pool_size,
895  max_element_size,
896  max_element_size);
897  S = Kokkos::View<scalar ***, Kokkos::LayoutRight, Kokkos::DefaultExecutionSpace>("A",
898  pool_size,
899  3 * max_element_size,
900  max_element_size);
901 
902  ElementMatrices<scalar, static_cast<int>(tbox::ElType::HEX8), 8> *elementMatrices_HEX8 =
903  new ElementMatrices<scalar, static_cast<int>(tbox::ElType::HEX8), 8>(domain, numPrimalDPN, pool_size);
904  ElementMatrices<scalar, static_cast<int>(tbox::ElType::TETRA4), 4> *elementMatrices_TETRA4 =
905  new ElementMatrices<scalar, static_cast<int>(tbox::ElType::TETRA4), 4>(domain, numPrimalDPN, pool_size);
906 
907  Kokkos::parallel_for(
908  numMyElems, KOKKOS_LAMBDA(const int i_elem) {
909  bool bContinue = true;
910 
911 #ifdef KOKKOS_ENABLE_DEPRECATED_CODE
912  int i_thread = Kokkos::DefaultExecutionSpace::thread_pool_rank();
913 #else
914  int i_thread = Kokkos::DefaultExecutionSpace::impl_thread_pool_rank();
915 #endif
916  const int e_size = domain->elementsList->getElementSize(i_elem);
917 
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));
921 
922  if (domain->elementsList->getElementType(i_elem) == static_cast<int>(tbox::ElType::HEX8))
923  {
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));
928  }
929 
930  else if (domain->elementsList->getElementType(i_elem) == static_cast<int>(tbox::ElType::TETRA4))
931  {
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));
936  }
937  else
938  bContinue = false;
939 
940  if (bContinue)
941  for (auto i = 0; i < e_size; ++i)
942  {
943  local_ordinal_type node_i = domain->elementsList->getElementNode(i_elem, i);
944  global_ordinal_type global_node_id_i = algebraic->map->mapNodesWO->getGlobalElement(node_i);
945  // Test if mapNodes ownes global_node_id_i
946  if (algebraic->map->mapNodes->isNodeGlobalElement(global_node_id_i))
947  {
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)
950  {
951  local_ordinal_type node_j = domain->elementsList->getElementNode(i_elem, j);
952 
953  if (numPrimalDPN == 4 || numPrimalDPN == 1)
954  {
955  local_ordinal_type cols = node_j;
956  scalar vals = L_e(i, j);
957 
958  local_L.sumIntoValues(local_node_id_i, &cols, 1, &vals, false, true);
959  }
960  if (numPrimalDPN == 4 || numPrimalDPN == 3)
961  {
962  for (auto k = 0; k < 3; ++k)
963  for (auto l = 0; l < 3; ++l)
964  {
965  local_ordinal_type cols = 3 * node_j + 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);
968  }
969  }
970  if (numPrimalDPN == 4)
971  {
972  for (auto k = 0; k < 3; ++k)
973  {
974  local_ordinal_type cols = node_j;
975  scalar vals = S_e(3 * i + k, j);
976  local_S.sumIntoValues(3 * local_node_id_i + k, &cols, 1, &vals, false, true);
977  }
978  }
979  }
980  }
981  }
982  });
983 
984  if (numPrimalDPN == 4 || numPrimalDPN == 1)
985  {
986  algebraic->matrices->L->fillComplete();
987  }
988  if (numPrimalDPN == 4 || numPrimalDPN == 3)
989  {
990  algebraic->matrices->K->fillComplete();
991  }
992  if (numPrimalDPN == 4)
993  {
994  algebraic->matrices->S->fillComplete();
995  }
996 
997  delete elementMatrices_HEX8;
998  delete elementMatrices_TETRA4;
999 
1000  if (numPrimalDPN == 4)
1001  {
1002  RCP<xpetra_bcrs_type> bA = rcp(new xpetra_bcrs_type(algebraic->map->thermomecBlockedMap, algebraic->map->thermomecBlockedMap, 8));
1003 
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));
1007 
1008  RCP<xpetra_wcrs_type> xwL = rcp(new xpetra_wcrs_type(xL));
1009  RCP<xpetra_wcrs_type> xwK = rcp(new xpetra_wcrs_type(xK));
1010  RCP<xpetra_wcrs_type> xwS = rcp(new xpetra_wcrs_type(xS));
1011 
1012  bA->setMatrix(0, 0, xwL);
1013  bA->setMatrix(1, 0, xwS);
1014  bA->setMatrix(1, 1, xwK);
1015 
1016  bA->fillComplete();
1017 
1018  algebraic->matrices->A = bA;
1019  }
1020  else if (numPrimalDPN == 3)
1021  {
1022  RCP<xpetra_crs_type> xK = rcp(new xpetra_tcrs_type(algebraic->matrices->K));
1023  algebraic->matrices->A = rcp(new xpetra_wcrs_type(xK));
1024  }
1025  else
1026  {
1027  RCP<xpetra_crs_type> xL = rcp(new xpetra_tcrs_type(algebraic->matrices->L));
1028  algebraic->matrices->A = rcp(new xpetra_wcrs_type(xL));
1029  }
1030 }
1031 
1040 template <typename scalar>
1042 {
1043 
1044  using Teuchos::RCP;
1045  using Teuchos::rcp;
1046  using Teuchos::tuple;
1047 
1048  typedef Kokkos::View<scalar ***, Kokkos::LayoutRight, Kokkos::DefaultExecutionSpace> elem_matrix_view_type;
1049 
1050  elem_matrix_view_type N, N_S;
1051 
1052  auto numMyElems = domain->elementsList->getElementNumber();
1053 
1054 #ifdef KOKKOS_ENABLE_DEPRECATED_CODE
1055  int pool_size = Kokkos::DefaultExecutionSpace::thread_pool_size();
1056 #else
1057  int pool_size = Kokkos::DefaultExecutionSpace::impl_thread_pool_size();
1058 #endif
1059  int max_element_size = 1;
1060  int max_element_size_S = 1;
1061 
1062  for (int i_elem = 0; i_elem < numMyElems; ++i_elem)
1063  {
1064  if (domain->elementsList->getElementType(i_elem) == static_cast<int>(tbox::ElType::HEX8))
1065  {
1066  max_element_size = 8;
1067  max_element_size_S = 4;
1068  break;
1069  }
1070  if (domain->elementsList->getElementType(i_elem) == static_cast<int>(tbox::ElType::TETRA4))
1071  {
1072  max_element_size = 4;
1073  max_element_size_S = 3;
1074  break;
1075  }
1076  }
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);
1079 
1080  ElementVectors<scalar, static_cast<int>(tbox::ElType::TRI3), 3> *elementVectors_TRI3 =
1081  new ElementVectors<scalar, static_cast<int>(tbox::ElType::TRI3), 3>(domain, numPrimalDPN, pool_size);
1082  ElementVectors<scalar, static_cast<int>(tbox::ElType::QUAD4), 4> *elementVectors_QUAD4 =
1083  new ElementVectors<scalar, static_cast<int>(tbox::ElType::QUAD4), 4>(domain, numPrimalDPN, pool_size);
1084  ElementVectors<scalar, static_cast<int>(tbox::ElType::TETRA4), 4> *elementVectors_TETRA4 =
1085  new ElementVectors<scalar, static_cast<int>(tbox::ElType::TETRA4), 4>(domain, numPrimalDPN, pool_size);
1086  ElementVectors<scalar, static_cast<int>(tbox::ElType::HEX8), 8> *elementVectors_HEX8 =
1087  new ElementVectors<scalar, static_cast<int>(tbox::ElType::HEX8), 8>(domain, numPrimalDPN, pool_size);
1088 
1089  RCP<tpetra_vector_type> bWO = rcp(new tpetra_vector_type(algebraic->map->mapDofsWO, true));
1090 
1091  auto bWO_view = bWO->getLocalViewHost();
1092 
1093  // ------------------------------------------
1094  //
1095  // volumetric heat sources
1096  //
1097  // ------------------------------------------
1098 
1099  auto numMySources = loads->sourcesList->getSourceNumber();
1100 
1101  for (auto i = 0; i < numMySources; ++i)
1102  {
1103  auto currentSourceSize = loads->sourcesList->getSourceSize(i);
1104  Kokkos::parallel_for(
1105  currentSourceSize, KOKKOS_LAMBDA(const int e1) {
1106  local_ordinal_type i_elem = loads->sourcesList->getSourceElement(i, e1);
1107 #ifdef KOKKOS_ENABLE_DEPRECATED_CODE
1108  int i_thread = Kokkos::DefaultExecutionSpace::thread_pool_rank();
1109 #else
1110  int i_thread = Kokkos::DefaultExecutionSpace::impl_thread_pool_rank();
1111 #endif
1112  const int e_size = domain->elementsList->getElementSize(i_elem);
1113 
1114  auto N_e = subview(N, i_thread, make_pair(0, e_size), make_pair(0, e_size));
1115 
1116  bool bContinue = true;
1117 
1118  if (domain->elementsList->getElementType(i_elem) == static_cast<int>(tbox::ElType::HEX8))
1119  {
1120  elementVectors_HEX8->compute(i_elem, i_thread);
1121  Kokkos::deep_copy(N_e, elementVectors_HEX8->getN(i_thread));
1122  }
1123  else if (domain->elementsList->getElementType(i_elem) == static_cast<int>(tbox::ElType::TETRA4))
1124  {
1125  elementVectors_TETRA4->compute(i_elem, i_thread);
1126  Kokkos::deep_copy(N_e, elementVectors_TETRA4->getN(i_thread));
1127  }
1128  else
1129  bContinue = false;
1130 
1131  if (bContinue)
1132  {
1133  for (auto ii = 0; ii < e_size; ++ii)
1134  {
1135  local_ordinal_type local_i = domain->elementsList->getElementNode(i_elem, ii);
1136  global_ordinal_type global_i = algebraic->map->mapNodesWO->getGlobalElement(local_i);
1137 
1138  // Test if mapDofs ownes global_i
1139  if (algebraic->map->mapNodes->isNodeGlobalElement(global_i))
1140  {
1141  scalar tmp = 0.;
1142 
1143  for (auto jj = 0; jj < e_size; ++jj)
1144  tmp += N_e(ii, jj);
1145 
1146  tmp *= loads->sourcesList->getSourceValue(i);
1147 
1148  Kokkos::atomic_fetch_add(&bWO_view(local_i, 0), tmp);
1149  }
1150  }
1151  }
1152  });
1153  }
1154 
1155  // ------------------------------------------
1156  //
1157  // von Neumann boundary conditions
1158  //
1159  // ------------------------------------------
1160 
1161  auto numMynBCs = loads->neumannList->getNeumannNumber();
1162 
1163  const size_t numNodes = algebraic->map->mapNodes->getGlobalNumElements();
1164  //const size_t numMyNodes = algebraic->map->mapNodes->getNodeNumElements();
1165  const size_t numMyNodesWO = algebraic->map->mapNodesWO->getNodeNumElements();
1166 
1167  for (auto i = 0; i < numMynBCs; ++i)
1168  {
1169  auto currentBCSize = loads->neumannList->getNeumannSize(i);
1170  Kokkos::parallel_for(
1171  currentBCSize, KOKKOS_LAMBDA(const int e1) {
1172  local_ordinal_type i_elem = loads->neumannList->getNeumannElement(i, e1);
1173 
1174 #ifdef KOKKOS_ENABLE_DEPRECATED_CODE
1175  int i_thread = Kokkos::DefaultExecutionSpace::thread_pool_rank();
1176 #else
1177  int i_thread = Kokkos::DefaultExecutionSpace::impl_thread_pool_rank();
1178 #endif
1179  const int e_size = domain->elementsList->getElementSize(i_elem);
1180 
1181  auto N_e = subview(N_S, i_thread, make_pair(0, e_size), make_pair(0, e_size));
1182 
1183  bool bContinue = true;
1184 
1185  if (domain->elementsList->getElementType(i_elem) == static_cast<int>(tbox::ElType::TRI3))
1186  {
1187  elementVectors_TRI3->compute(i_elem, i_thread);
1188  Kokkos::deep_copy(N_e, elementVectors_TRI3->getN(i_thread));
1189  }
1190  else if (domain->elementsList->getElementType(i_elem) == static_cast<int>(tbox::ElType::QUAD4))
1191  {
1192  elementVectors_QUAD4->compute(i_elem, i_thread);
1193  Kokkos::deep_copy(N_e, elementVectors_QUAD4->getN(i_thread));
1194  }
1195  else
1196  bContinue = false;
1197 
1198  if (bContinue)
1199  {
1200  for (auto j = 0; j < numPrimalDPN; ++j)
1201  {
1202  scalar which_dof = loads->neumannList->getNeumannDof(i, j);
1203  scalar value = loads->neumannList->getNeumannValue(i, j);
1204  if (which_dof == (scalar)1)
1205  {
1206  for (auto ii = 0; ii < e_size; ++ii)
1207  {
1208  local_ordinal_type local_i;
1209  if (j == 3)
1210  local_i = domain->elementsList->getElementNode(i_elem, ii);
1211  else
1212  {
1213  if (numPrimalDPN == 4)
1214  local_i = numMyNodesWO + 3 * domain->elementsList->getElementNode(i_elem, ii) + j;
1215  else
1216  local_i = 3 * domain->elementsList->getElementNode(i_elem, ii) + j;
1217  }
1218 
1219  scalar tmp = 0.;
1220 
1221  for (auto jj = 0; jj < e_size; ++jj)
1222  tmp += N_e(ii, jj);
1223 
1224  tmp *= value;
1225 
1226  Kokkos::atomic_fetch_add(&bWO_view(local_i, 0), tmp);
1227  }
1228  }
1229  }
1230  }
1231  });
1232  }
1233 
1234  // ------------------------------------------
1235  //
1236  // Dirichlet boundary conditions
1237  //
1238  // ------------------------------------------
1239 
1240  Tpetra::Export<> exportBWO(algebraic->map->mapDofsWO, algebraic->map->mapDofs);
1241  algebraic->vectors->b->doExport(*bWO, exportBWO, Tpetra::ADD);
1242 
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();
1249 
1250  auto numMyBCs = loads->dirichletList->getDirichletNumber();
1251  for (auto i = 0; i < numMyBCs; ++i)
1252  {
1253  auto currentBCSize = loads->dirichletList->getDirichletSize(i);
1254  for (int i_node = 0; i_node < currentBCSize; ++i_node)
1255  {
1256  global_ordinal_type n1 = loads->dirichletList->getDirichletNode(i, i_node);
1257 
1258  if (algebraic->map->mapNodes->isNodeGlobalElement(n1))
1259  for (auto j = 0; j < numPrimalDPN; ++j)
1260  {
1261  scalar which_dof = loads->dirichletList->getDirichletDof(i, j);
1262 
1263  if (which_dof == (scalar)1)
1264  {
1265  global_ordinal_type global_index, global_index_0, global_index_r0;
1266 
1267  bool thermo = false;
1268 
1269  if (j == 3)
1270  {
1271  global_index_0 = 0;
1272  global_index_r0 = n1;
1273 
1274  thermo = true;
1275  }
1276  else
1277  {
1278  if (numPrimalDPN == 4)
1279  global_index_0 = numNodes;
1280  else
1281  global_index_0 = 0;
1282  global_index_r0 = 3 * n1 + j;
1283  }
1284 
1285  global_index = global_index_0 + global_index_r0;
1286 
1287  Teuchos::Array<global_ordinal_type> global_indices;
1288 
1289  if (thermo)
1290  {
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);
1294 
1295  // Replace the initialized elements by zero
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)));
1300 
1301  // Put one on the diagonal element
1302  algebraic->matrices->L->replaceGlobalValues(global_index,
1303  tuple<global_ordinal_type>(global_index),
1304  tuple<scalar>(scalar(1.0)));
1305  }
1306  else
1307  {
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);
1311 
1312  // Replace the initialized elements by zero
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)));
1317 
1318  // Put one on the diagonal element
1319  algebraic->matrices->K->replaceGlobalValues(global_index_r0,
1320  tuple<global_ordinal_type>(global_index_r0),
1321  tuple<scalar>(scalar(1.0)));
1322 
1323  if (numPrimalDPN == 4)
1324  {
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);
1328 
1329  // Replace the initialized elements by zero
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)));
1334  }
1335  }
1336 
1337  // Put the corresponding value in the right hand side
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))));
1342  }
1343  }
1344  }
1345  }
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();
1352 
1353  /*
1354  if (scaled)
1355  {
1356  tpetra_vector_type diag(algebraic->matrices->A->getRowMap());
1357  algebraic->matrices->A->getLocalDiagCopy(diag);
1358 
1359  Teuchos::ArrayRCP<scalar> diag_entries;
1360  diag_entries = diag.getDataNonConst();
1361  for (size_t i = 0; i < diag_entries.size(); ++i)
1362  diag_entries[i] = 1. / diag_entries[i];
1363 
1364  Teuchos::ArrayRCP<scalar> rhs_entries;
1365  rhs_entries = algebraic->vectors->b->getDataNonConst();
1366  for (size_t i = 0; i < rhs_entries.size(); ++i)
1367  rhs_entries[i] = rhs_entries[i] * diag_entries[i];
1368 
1369  algebraic->matrices->A->leftScale(diag);
1370 
1371  if (domain->contactsList->getContactNumber() >= 1)
1372  algebraic->matrices->B_T->leftScale(diag);
1373  }
1374  */
1375  // ------------------------------------------
1376  //
1377  // weights computation
1378  //
1379  // ------------------------------------------
1380 
1381  auto numMyWeightRegions = loads->weightsList->getWeightRegionsNumber();
1382 
1383  for (auto i = 0; i < numMyWeightRegions; ++i)
1384  {
1385  auto currentWeightRegionSize = loads->weightsList->getWeightRegionSize(i);
1386  Kokkos::parallel_for(
1387  currentWeightRegionSize, KOKKOS_LAMBDA(const int i_node) {
1388  global_ordinal_type n1 = loads->weightsList->getNode(i, i_node);
1389 
1390  if (algebraic->map->mapNodes->isNodeGlobalElement(n1))
1391  {
1392  for (auto j = 0; j < numPrimalDPN; ++j)
1393  {
1394  if (loads->weightsList->getWeightDof(i, j))
1395  {
1396  global_ordinal_type global_index;
1397 
1398  if (j == 3)
1399  {
1400  global_index = n1;
1401  }
1402  else
1403  {
1404  if (numPrimalDPN == 4)
1405  global_index = numNodes + 3 * n1 + j;
1406  else
1407  global_index = 3 * n1 + j;
1408  }
1409 
1410  algebraic->vectors->weights->sumIntoGlobalValue(global_index,
1411  (scalar(loads->weightsList->getWeightValue(i, j))),
1412  true);
1413  }
1414  }
1415  }
1416  });
1417  }
1418 
1419  delete elementVectors_TRI3;
1420  delete elementVectors_QUAD4;
1421  delete elementVectors_TETRA4;
1422  delete elementVectors_HEX8;
1423 }
1424 }; // namespace katoptron
katoptron::Domain
Class which is used to store all the information related to the discretized domain:
Definition: Domain.h:26
katoptron::DiscreteProblem::xpetra_bmap_type
Map::xpetra_bmap_type xpetra_bmap_type
Definition: DiscreteProblem.h:58
katoptron::DiscreteProblem::xpetra_tcrs_type
Matrices< scalar >::xpetra_tcrs_type xpetra_tcrs_type
Definition: DiscreteProblem.h:70
katoptron::DiscreteProblem::xpetra_wcrs_type
Matrices< scalar >::xpetra_wcrs_type xpetra_wcrs_type
Definition: DiscreteProblem.h:68
ElementMatrices
Class used to compute the element matrices:
Definition: ElementMatrices.h:46
katoptron::DiscreteProblem::global_ordinal_type
Map::global_ordinal_type global_ordinal_type
Definition: DiscreteProblem.h:51
katoptron::DiscreteProblem::xpetra_tmvector_type
Vectors< scalar >::xpetra_tmvector_type xpetra_tmvector_type
Definition: DiscreteProblem.h:64
ElementMatrices::L
view_type_3D L
Definition: ElementMatrices.h:52
katoptron::DiscreteProblem::xpetra_bmvector_type
Vectors< scalar >::xpetra_bmvector_type xpetra_bmvector_type
Definition: DiscreteProblem.h:62
katoptron::DiscreteProblem::xpetra_tmap_type
Map::xpetra_tmap_type xpetra_tmap_type
Definition: DiscreteProblem.h:57
katoptron::DiscreteProblem::local_ordinal_type
Map::local_ordinal_type local_ordinal_type
Definition: DiscreteProblem.h:50
katoptron::DiscreteProblem::export_type
Tpetra::Export< local_ordinal_type, global_ordinal_type > export_type
Definition: DiscreteProblem.h:53
katoptron::DiscreteProblem::tpetra_map_type
Map::tpetra_map_type tpetra_map_type
Definition: DiscreteProblem.h:55
katoptron::Algebraic
Class which is used to store Teuchos::RCP to the algebraic information of the problem:
Definition: Algebraic.h:31
katoptron
katoptron namespace
Definition: Algebraic.h:18
katoptron::DiscreteProblem::computeMatrices
void computeMatrices(void)
Compute the primal matrix.
Definition: DiscreteProblem.hpp:831
katoptron::DiscreteProblem::computeLoads
void computeLoads(bool scaled=false)
Compute the rhs.
Definition: DiscreteProblem.hpp:1041
katoptron::DiscreteProblem::computeContactMatrices
void computeContactMatrices(MPI_Comm comm)
Compute the contacts matrices.
Definition: DiscreteProblem.hpp:818
katoptron::DiscreteProblem::tpetra_vector_type
Vectors< scalar >::tpetra_vector_type tpetra_vector_type
Definition: DiscreteProblem.h:60
katoptron::LinearSolver::getTimers
static fwk::Timers & getTimers()
Manage timers.
Definition: LinearSolver.h:119
katoptron::Loads
Class which includes all the loads, boundary conditions, preloads, and weights.
Definition: Loads.h:21
katoptron::DiscreteProblem::xpetra_bcrs_type
Matrices< scalar >::xpetra_bcrs_type xpetra_bcrs_type
Definition: DiscreteProblem.h:69
katoptron::Problem
Class which is used to specify in Python the thermomechanical to solve.
Definition: wProblem.h:19
katoptron::Map
Class which includes all the Trilinos maps (Tpetra maps and Xpetra maps) used in the simulation.
Definition: Map.h:21
katoptron::DiscreteProblem::createBlockMVector
void createBlockMVector(bool write)
Create the block vectors.
Definition: DiscreteProblem.hpp:714
katoptron::DiscreteProblem::DiscreteProblem
DiscreteProblem(Problem &pbl, size_t _numPrimalDPN, Teuchos::RCP< Teuchos::ParameterList > randomParams, Kokkos::View< scalar *, Kokkos::LayoutLeft > m_rv_values)
DiscreteProblem constructor.
Definition: DiscreteProblem.hpp:28
katoptron::DiscreteProblem::updateActivity
Mask< scalar > updateActivity(int active_set_iteration)
Update the activity of the Lagrange multipliers and return a Mask which stores true for the converged...
Definition: DiscreteProblem.hpp:56
ElementVectors
Class used to compute the matrix associated to the integration of load on surface .
Definition: ElementVectors.h:27
katoptron::DiscreteProblem::createBlockMatrix
void createBlockMatrix(scalar gamma, bool write, bool merge)
Create the block matrix.
Definition: DiscreteProblem.hpp:479