/[escript]/trunk/finley/src/Mesh_getTrilinosGraph.cpp
ViewVC logotype

Contents of /trunk/finley/src/Mesh_getTrilinosGraph.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6939 - (show annotations)
Mon Jan 20 03:37:18 2020 UTC (4 months, 2 weeks ago) by uqaeller
File size: 3671 byte(s)
Updated the copyright header.


1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2020 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Apache License, version 2.0
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014-2017 by Centre for Geoscience Computing (GeoComp)
14 * Development from 2019 by School of Earth and Environmental Sciences
15 **
16 *****************************************************************************/
17
18 #ifdef ESYS_HAVE_TRILINOS
19 #include "FinleyDomain.h"
20 #include "IndexList.h"
21
22 #include <boost/scoped_array.hpp>
23
24 using namespace esys_trilinos;
25
26 namespace finley {
27
28 const_TrilinosGraph_ptr FinleyDomain::getTrilinosGraph(bool reducedOrder) const
29 {
30 const_TrilinosGraph_ptr out;
31 // make sure that the requested graph is available
32 if (reducedOrder) {
33 if (m_reducedGraph.is_null())
34 m_reducedGraph.reset(createTrilinosGraph(reducedOrder));
35 out = m_reducedGraph;
36 } else {
37 if (m_fullGraph.is_null())
38 m_fullGraph.reset(createTrilinosGraph(reducedOrder));
39 out = m_fullGraph;
40 }
41 return out;
42 }
43
44 GraphType* FinleyDomain::createTrilinosGraph(bool reducedOrder) const
45 {
46 index_t myNumTargets;
47 index_t numTargets;
48 const index_t* target;
49 const_TrilinosMap_ptr rowMap;
50 const_TrilinosMap_ptr colMap;
51 if (reducedOrder) {
52 myNumTargets = m_nodes->getNumReducedDegreesOfFreedom();
53 numTargets = m_nodes->getNumReducedDegreesOfFreedomTargets();
54 target = m_nodes->borrowTargetReducedDegreesOfFreedom();
55 rowMap = m_nodes->trilinosReducedRowMap;
56 colMap = m_nodes->trilinosReducedColMap;
57 } else {
58 myNumTargets = m_nodes->getNumDegreesOfFreedom();
59 numTargets = m_nodes->getNumDegreesOfFreedomTargets();
60 target = m_nodes->borrowTargetDegreesOfFreedom();
61 rowMap = m_nodes->trilinosRowMap;
62 colMap = m_nodes->trilinosColMap;
63 }
64
65 boost::scoped_array<IndexList> indexList(new IndexList[numTargets]);
66
67 #pragma omp parallel
68 {
69 // insert contributions from element matrices into columns in
70 // index list
71 IndexList_insertElements(indexList.get(), m_elements, reducedOrder,
72 target, reducedOrder, target);
73 IndexList_insertElements(indexList.get(), m_faceElements,
74 reducedOrder, target, reducedOrder, target);
75 IndexList_insertElements(indexList.get(), m_contactElements,
76 reducedOrder, target, reducedOrder, target);
77 IndexList_insertElements(indexList.get(), m_points, reducedOrder,
78 target, reducedOrder, target);
79 }
80
81 Teuchos::ArrayRCP<size_t> rowPtr(myNumTargets + 1);
82 for (size_t i = 0; i < myNumTargets; i++) {
83 rowPtr[i+1] = rowPtr[i] + indexList[i].count(0, numTargets);
84 }
85
86 Teuchos::ArrayRCP<LO> colInd(rowPtr[myNumTargets]);
87
88 #pragma omp parallel for
89 for (index_t i = 0; i < myNumTargets; i++) {
90 indexList[i].toArray(&colInd[rowPtr[i]], 0, numTargets, 0);
91 std::sort(&colInd[rowPtr[i]], &colInd[rowPtr[i+1]]);
92 }
93
94 GraphType* graph = new GraphType(rowMap, colMap, rowPtr, colInd);
95 Teuchos::RCP<Teuchos::ParameterList> params = Teuchos::parameterList();
96 params->set("Optimize Storage", true);
97 graph->fillComplete(rowMap, rowMap, params);
98 return graph;
99 }
100
101 } // namespace finley
102
103 #endif // ESYS_HAVE_TRILINOS
104

  ViewVC Help
Powered by ViewVC 1.1.26