FromVertexMatching.cpp 9.87 KB
Newer Older
1
2
//! \file 
//
3
4
5
6
7
8
9
10
//
//  FromVertexMatching.cpp
//  HappyHeart
//
//  Created by Sebastien Gilles on 14/12/15.
//  Copyright © 2015 Inria. All rights reserved.
//

11
#include "HappyHeart/ThirdParty/Wrappers/Petsc/Matrix/MatrixPattern.hpp"
12

13
#include "HappyHeart/Core/NumberingSubset/NumberingSubset.hpp"
14

15
16
17
18
#include "HappyHeart/FiniteElement/FiniteElementSpace/FEltSpace.hpp"
#include "HappyHeart/FiniteElement/FiniteElementSpace/Internal/DofProgramWiseIndexListPerVertexCoordIndex.hpp"
#include "HappyHeart/FiniteElement/FiniteElementSpace/Internal/DofProgramWiseIndexListPerVertexCoordIndexManager.hpp"
#include "HappyHeart/FiniteElement/FiniteElementSpace/GodOfDof.hpp"
19

20
#include "HappyHeart/OperatorInstances/NonConformInterpolator/FromVertexMatching.hpp"
21
22
23
24
25
26
27
28
29
30


namespace HappyHeart
{
    
    
    namespace NonConformInterpolatorNS
    {
        
        
31
        void FromVertexMatching::Construct(const MeshNS::InterpolationNS::VertexMatching& vertex_matching,
32
33
                                           const unsigned int source_index,
                                           const unsigned int target_index)
34
        {
35
36
37
38
39
40
41
42
43
            auto& init_vertex_matching_manager =
                Internal::FEltSpaceNS::DofProgramWiseIndexListPerVertexCoordIndexManager::GetInstance();

            decltype(auto) source =
                init_vertex_matching_manager.GetDofProgramWiseIndexListPerVertexCoordIndex(source_index);
            
            decltype(auto) target =
                init_vertex_matching_manager.GetDofProgramWiseIndexListPerVertexCoordIndex(target_index);
            
44
            interpolation_matrix_ = std::make_unique<GlobalMatrix>(target.GetNumberingSubset(),
45
                                                                   source.GetNumberingSubset());
46
            
47
            decltype(auto) source_dof_index_list_per_coord_index = source.GetDofProgramWiseIndexPerCoordIndex();
48
            decltype(auto) target_dof_list_per_proc_wise_coord_index = target.ComputeProcessorWiseStorage();
49
50
51
            
            const auto Nprocessor_wise_source_dof = source.ComputeNprocessorWisedof();
            const auto Nprocessor_wise_target_dof = target.ComputeNprocessorWisedof();
52
            
53
54
            const auto Nprogram_wise_source_dof = source.NprogramWisedof();
            const auto Nprogram_wise_target_dof = target.NprogramWisedof();
55
            
56
57
            assert(Nprogram_wise_source_dof == Nprogram_wise_target_dof);
            assert(!source_dof_index_list_per_coord_index.empty());
58
            
59
            std::vector<std::vector<PetscInt>> non_zero_pattern_per_line(Nprocessor_wise_target_dof);
60
61
62
63
64
            
            const auto dumb_value = static_cast<PetscInt>(-1);
            
            std::vector<PetscInt> target_processor_to_program_wise(Nprocessor_wise_target_dof,
                                                                   dumb_value);
65
            
66
            decltype(auto) target_numbering_subset = target.GetNumberingSubset();
67

68
69
70
71
            const auto complete_dof_list_target = target.GetFEltSpace().GetProcessorWiseDofList(target_numbering_subset);
            
            const auto end_dof_list_target = complete_dof_list_target.cend();
            
72
            for (const auto& pair : target_dof_list_per_proc_wise_coord_index)
73
74
75
            {
                // Use geometric data to make source and target vertices match.
                // Consider only Coords handled processor-wisely.
76
77
                const auto target_coords_vertex_index = pair.first;
                const auto source_coords_vertex_index = vertex_matching.FindSourceIndex(target_coords_vertex_index);
78
79
                
                // Then read the associated dof data.
80
                const auto it_source = source_dof_index_list_per_coord_index.find(source_coords_vertex_index);
81
82
83
                assert(it_source != source_dof_index_list_per_coord_index.cend());
                    
                const auto& source_dof_list = it_source->second;
84
                const auto& target_dof_list = pair.second;
85
86
87
88
89
90
                
                const auto Ndof_for_vertex = source_dof_list.size();
                assert(Ndof_for_vertex == target_dof_list.size());
                
                for (auto j = 0ul; j < Ndof_for_vertex; ++j)
                {
91
                    // \todo #775 Ordering of the unknowns shouldn't be implicitly assumed.
92
93
                    const auto target_program_wise_index = target_dof_list[j];
                    const auto source_program_wise_index = source_dof_list[j];
94
                    
95
96
                    // Identify the dof objects: we need respectively processor- and program-wise indexes for row
                    // and columns.
97
                    const auto it_target =
98
99
                        std::find_if(complete_dof_list_target.cbegin(),
                                     end_dof_list_target,
100
                                     [target_program_wise_index, &target_numbering_subset](const auto& dof_ptr)
101
102
                                     {
                                         assert(!(!dof_ptr));
103
                                         return dof_ptr->GetProgramWiseIndex(target_numbering_subset) == target_program_wise_index;
104
105
106
107
                                     });
                    
                    // This might happen: a Coords may exist processor-wise but bear no dofs (they are in this case born
                    // by another processor).
108
                    if (it_target == end_dof_list_target)
109
                        continue;
110
111
112
                    
                    assert(!(!*it_target));
                    const auto& target_dof = *(*it_target);
113

114
                    const auto target_processor_wise_index = target_dof.GetProcessorWiseOrGhostIndex(target_numbering_subset);
115
                    
116
117
                    assert(target_processor_wise_index < Nprocessor_wise_target_dof);
                    assert(source_program_wise_index < Nprogram_wise_source_dof);
118
                    
119
                    // I reconstitute here a processor-2-program-wise array, usually handled through Dof object
120
                    // (can't be used here - see Internal::DofProgramWiseIndexListPerVertexCoordIndex for more about it).
121
122
123
124
125
                    non_zero_pattern_per_line[static_cast<std::size_t>(target_processor_wise_index)] =
                        { static_cast<PetscInt>(source_program_wise_index) };
                    
                    target_processor_to_program_wise[static_cast<std::size_t>(target_processor_wise_index)] =
                        static_cast<PetscInt>(target_program_wise_index);
126
127
                }
            }
128
            
129
130
131
132
133
134
135
136
            assert(std::none_of(target_processor_to_program_wise.cbegin(),
                                target_processor_to_program_wise.cend(),
                                [dumb_value](auto value)
                                {
                                    return value == dumb_value;
                                }));
            
            
137
138
139
140
            matrix_pattern_ = std::make_unique<Wrappers::Petsc::MatrixPattern>(non_zero_pattern_per_line);
            
            decltype(auto) matrix_pattern = GetMatrixPattern();
            
141
            auto& interpolation_matrix = *interpolation_matrix_;
142
            
143
            const auto& mpi = source.GetFEltSpace().MpiHappyHeart();
144
145
146
            
            if (mpi.Nprocessor<int>() == 1)
            {
GILLES Sebastien's avatar
GILLES Sebastien committed
147
148
149
150
                assert(Nprocessor_wise_source_dof == Nprocessor_wise_target_dof);
                
                interpolation_matrix.InitSequentialMatrix(static_cast<unsigned int>(Nprocessor_wise_target_dof),
                                                          static_cast<unsigned int>(Nprocessor_wise_source_dof),
151
                                                          matrix_pattern,
152
                                                          mpi,
153
154
155
156
                                                          __FILE__, __LINE__);
            }
            else
            {
GILLES Sebastien's avatar
GILLES Sebastien committed
157
                assert(Nprogram_wise_source_dof == Nprogram_wise_target_dof);
158
                
GILLES Sebastien's avatar
GILLES Sebastien committed
159
160
                interpolation_matrix.InitParallelMatrix(static_cast<unsigned int>(Nprocessor_wise_target_dof),
                                                        static_cast<unsigned int>(Nprocessor_wise_source_dof),
161
162
163
164
165
166
167
                                                        Nprogram_wise_target_dof,
                                                        Nprogram_wise_source_dof,
                                                        matrix_pattern,
                                                        mpi,
                                                        __FILE__, __LINE__);
            }

168
            
169
            // Fill the interpolation matrix.
170
            std::vector<PetscScalar> one { 1. };
171
172
173
174
175
176
177
178
                    
            assert(non_zero_pattern_per_line.size() == Nprocessor_wise_target_dof);
            assert(target_processor_to_program_wise.size() == Nprocessor_wise_target_dof);
            
            std::vector<PetscInt> source_index_list;
            source_index_list.reserve(Nprocessor_wise_target_dof);
            
            for (decltype(auto) item : non_zero_pattern_per_line)
179
            {
180
181
182
183
184
                assert(item.size() == 1ul);
                source_index_list.push_back(item.back());
            }
            
            std::vector<PetscScalar> one_list(Nprocessor_wise_target_dof, 1.);
185

186
187
188
189
190
191
192
            for (auto i = 0ul; i < Nprocessor_wise_target_dof; ++i)
            {
                interpolation_matrix.SetValue(target_processor_to_program_wise[i],
                                              source_index_list[i],
                                              1.,
                                              INSERT_VALUES,
                                              __FILE__, __LINE__);
193
194
195
            }
            
            interpolation_matrix.Assembly(__FILE__, __LINE__);
196
197
198
        }
        
        
199
200
201
202
    } // namespace NonConformInterpolatorNS


} // namespace HappyHeart