1 #ifndef KATOPTRON_MORTAR_H
2 #define KATOPTRON_MORTAR_H
4 #include "Tpetra_Vector.hpp"
5 #include "Kokkos_ViewFactory.hpp"
6 #include "Tpetra_Map.hpp"
8 #include "Tpetra_CrsGraph.hpp"
9 #include <Epetra_Map.h>
11 #include "Tpetra_ConfigDefs.hpp"
13 #include "Tpetra_Map.hpp"
14 #include "Tpetra_MultiVector.hpp"
15 #include "Tpetra_Vector.hpp"
16 #include "Tpetra_CrsGraph.hpp"
17 #include "Tpetra_CrsMatrix.hpp"
25 Teuchos::RCP<Tpetra::CrsMatrix<>>
compute_B(Teuchos::RCP<katoptron::ElementsList> elementsList,
26 Teuchos::RCP<katoptron::NodesList> nodesList,
27 Teuchos::RCP<katoptron::ContactsList> contactsList,
28 Teuchos::RCP<katoptron::Map> map, MPI_Comm comm,
bool SignoriniCheck);
30 template <
typename scalar>
32 Teuchos::RCP<katoptron::ElementsList> elementsList,
33 Teuchos::RCP<katoptron::NodesList> nodesList,
34 Teuchos::RCP<katoptron::ContactsList> contactsList,
35 Teuchos::RCP<katoptron::Map> map, MPI_Comm comm)
37 Teuchos::RCP<Tpetra::CrsMatrix<>> B_double =
compute_B(elementsList, nodesList, contactsList, map, comm,
true);
38 Teuchos::RCP<Tpetra::CrsMatrix<>> B_G_double;
39 if (contactsList->isSignorini())
40 B_G_double =
compute_B(elementsList, nodesList, contactsList, map, comm,
false);
42 B_G_double = B_double;
47 size_t maxNumEntPerRow = B_G_double->getNodeMaxNumRowEntries();
48 size_t maxNumEntPerRowT = 200;
50 matrices->B = rcp(
new Tpetra::CrsMatrix<scalar, local_ordinal_type, global_ordinal_type>(map->mapLagrangeDofs, maxNumEntPerRow));
51 matrices->B_2 = rcp(
new Tpetra::CrsMatrix<scalar, local_ordinal_type, global_ordinal_type>(map->mapLagrangeDofs, maxNumEntPerRow));
52 matrices->B_G = rcp(
new Tpetra::CrsMatrix<scalar, local_ordinal_type, global_ordinal_type>(map->mapLagrangeDofs, maxNumEntPerRow));
53 matrices->B_T = rcp(
new Tpetra::CrsMatrix<scalar, local_ordinal_type, global_ordinal_type>(map->mapDofs, maxNumEntPerRowT));
54 matrices->C = rcp(
new Tpetra::CrsMatrix<scalar, local_ordinal_type, global_ordinal_type>(map->mapLagrangeDofs, 1));
55 matrices->Cb = rcp(
new Tpetra::CrsMatrix<scalar, local_ordinal_type, global_ordinal_type>(map->mapLagrangeDofs, 1));
57 for (
auto i = 0; i < map->mapLagrangeDofs->getNodeNumElements(); ++i)
59 if (map->mapLagrangeDofs->isNodeLocalElement(i))
64 Teuchos::ArrayView<const local_ordinal_type> local_indices;
65 Teuchos::ArrayView<const double> values;
67 B_double->getLocalRowView(local_index, local_indices, values);
69 for (
auto k = 0; k < local_indices.size(); ++k)
71 global_ordinal_type column_index = B_double->getColMap()->getGlobalElement(local_indices[k]);
73 matrices->B->insertGlobalValues(global_index,
74 Teuchos::tuple<global_ordinal_type>(column_index),
75 Teuchos::tuple<scalar>((scalar)values[k]));
77 matrices->B_2->insertGlobalValues(global_index,
78 Teuchos::tuple<global_ordinal_type>(column_index),
79 Teuchos::tuple<scalar>((scalar)values[k]));
81 matrices->B_T->insertGlobalValues(column_index,
82 Teuchos::tuple<global_ordinal_type>(global_index),
83 Teuchos::tuple<scalar>((scalar)values[k]));
86 B_G_double->getLocalRowView(local_index, local_indices, values);
88 for (
auto k = 0; k < local_indices.size(); ++k)
90 global_ordinal_type column_index = B_G_double->getColMap()->getGlobalElement(local_indices[k]);
92 matrices->B_G->insertGlobalValues(global_index,
93 Teuchos::tuple<global_ordinal_type>(column_index),
94 Teuchos::tuple<scalar>((scalar)values[k]));
97 matrices->C->insertGlobalValues(global_index, Teuchos::tuple<global_ordinal_type>(global_index), Teuchos::tuple<scalar>((scalar)0.));
98 matrices->Cb->insertGlobalValues(global_index, Teuchos::tuple<global_ordinal_type>(global_index), Teuchos::tuple<scalar>((scalar)1.));
102 matrices->B->fillComplete(map->mapDofs, map->mapLagrangeDofs);
103 matrices->B_2->fillComplete(map->mapDofs, map->mapLagrangeDofs);
104 matrices->B_G->fillComplete(map->mapDofs, map->mapLagrangeDofs);
105 matrices->B_T->fillComplete(map->mapLagrangeDofs, map->mapDofs);
106 matrices->C->fillComplete();
107 matrices->Cb->fillComplete();
109 #endif //KATOPTRON_MORTAR_H