Skip to content

Commit ad1ee5e

Browse files
mikaemmikaem
committed
Compatible with fenics 2019.1
1 parent 99930d8 commit ad1ee5e

File tree

4 files changed

+105
-105
lines changed

4 files changed

+105
-105
lines changed

‎fenicstools/fem/cr_divergence.cpp‎

Lines changed: 27 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,7 @@ cfg['library_dirs'] = dolfin_pc['library_dirs']
2424

2525
using namespace dolfin;
2626

27-
// Base case for all divergence computations.
27+
// Base case for all divergence computations.
2828
// Compute divergence of vector field u.
2929
void cr_divergence(Function& divu, const Function& u)
3030
{
@@ -36,7 +36,7 @@ void cr_divergence(Function& divu, const Function& u)
3636
std::shared_ptr<const GenericVector> U = u.vector();
3737
std::shared_ptr<GenericVector> DIVU = divu.vector();
3838

39-
// Figure out about the local dofs of DG0
39+
// Figure out about the local dofs of DG0
4040
std::pair<std::size_t, std::size_t>
4141
first_last_dof = DG0_dofmap->ownership_range();
4242
std::size_t first_dof = first_last_dof.first;
@@ -48,7 +48,7 @@ void cr_divergence(Function& divu, const Function& u)
4848

4949
// Get topological dimension so that we know what Facet is
5050
const Mesh mesh = *u.function_space()->mesh();
51-
std::size_t tdim = mesh.topology().dim();
51+
std::size_t tdim = mesh.topology().dim();
5252
// Get the info on length of u and gdim for the dot product
5353
std::size_t gdim = mesh.geometry().dim();
5454
std::size_t udim = u.value_dimension(0); // Rank 1!
@@ -65,10 +65,10 @@ void cr_divergence(Function& divu, const Function& u)
6565
{
6666
Point cell_mp = cell->midpoint();
6767
double cell_volume = cell->volume();
68-
68+
6969
// Dofs of CR on all facets of the cell, global order
7070
auto facets_dofs = CR1_dofmap->cell_dofs(cell->index());
71-
71+
7272
double cell_integral = 0;
7373
std::size_t local_facet_index = 0;
7474
for(FacetIterator facet(*cell); !facet.end(); ++facet)
@@ -79,7 +79,7 @@ void cr_divergence(Function& divu, const Function& u)
7979
else if(tdim == 3)
8080
facet_measure = Face(mesh, facet->index()).area();
8181
// Tdim 1 will not happen because CR is not defined there
82-
82+
8383
Point facet_normal = facet->normal();
8484

8585
// Flip the normal if it is not outer already
@@ -115,7 +115,7 @@ void cr_divergence(Function& divu, const Function& u)
115115
void cr_divergence2(Function& divu, const Function& u,
116116
const FunctionSpace& DGscalar, const FunctionSpace& CRvector)
117117
{
118-
// Divu scalar from u vector
118+
// Divu scalar from u vector
119119
std::size_t u_rank = u.value_rank();
120120
std::size_t divu_rank = divu.value_rank();
121121
if((divu_rank == 0) and (u_rank == 1))
@@ -139,10 +139,10 @@ void cr_divergence2(Function& divu, const Function& u,
139139
std::shared_ptr<const GenericVector> U = u.vector();
140140
std::shared_ptr<GenericVector> DIVU = divu.vector();
141141

142-
// Get dofmaps
142+
// Get dofmaps
143143
std::shared_ptr<const FunctionSpace> DGvector = divu.function_space();
144144
std::shared_ptr<const FunctionSpace> CRtensor = u.function_space();
145-
145+
146146
std::size_t len_divu = divu.value_dimension(0);
147147
// With scalar U can be extracted only once
148148
std::vector<double> U_values;
@@ -163,7 +163,7 @@ void cr_divergence2(Function& divu, const Function& u,
163163
CRi_rows = CRvector[i]->dofmap()->dofs();
164164
std::size_t m = CRi_rows.size();
165165
U_i->zero();
166-
U_i->set(U_values.data(), m, CRi_rows.data());
166+
U_i->set(U_values.data(), m, CRi_rows.data());
167167
U_i->apply("insert");
168168
}
169169
// For tensor, U_i represents T_ij, (T_i0, T_i1, T_i2)
@@ -183,7 +183,7 @@ void cr_divergence2(Function& divu, const Function& u,
183183
std::vector<dolfin::la_index>
184184
CRj_rows = CRvector[j]->dofmap()->dofs();
185185

186-
U_i->set(U_values.data(), m, CRj_rows.data());
186+
U_i->set(U_values.data(), m, CRj_rows.data());
187187
U_i->apply("insert");
188188
}
189189
}
@@ -197,22 +197,22 @@ void cr_divergence2(Function& divu, const Function& u,
197197
std::size_t m = DGi_rows.size();
198198
std::vector<double> DIVU_i_values;
199199
DIVU_i->get_local(DIVU_i_values);
200-
DIVU->set(DIVU_i_values.data(), m, DGi_rows.data());
200+
DIVU->set(DIVU_i_values.data(), m, DGi_rows.data());
201201
DIVU->apply("insert");
202202
}
203203
}
204204
}
205-
// Base case for all divergence computations.
205+
// Base case for all divergence computations.
206206
// Compute divergence of vector field u.
207207
void cr_divergence_matrix(GenericMatrix& M, GenericMatrix& A,
208-
const FunctionSpace& DGscalar,
208+
const FunctionSpace& DGscalar,
209209
const FunctionSpace& CRvector)
210210
{
211211
std::shared_ptr<const GenericDofMap>
212212
CR1_dofmap = CRvector.dofmap(),
213213
DG0_dofmap = DGscalar.dofmap();
214214

215-
// Figure out about the local dofs of DG0
215+
// Figure out about the local dofs of DG0
216216
std::pair<std::size_t, std::size_t>
217217
first_last_dof = DG0_dofmap->ownership_range();
218218
std::size_t first_dof = first_last_dof.first;
@@ -221,25 +221,25 @@ void cr_divergence_matrix(GenericMatrix& M, GenericMatrix& A,
221221

222222
// Get topological dimension so that we know what Facet is
223223
const Mesh mesh = *DGscalar.mesh();
224-
std::size_t tdim = mesh.topology().dim();
224+
std::size_t tdim = mesh.topology().dim();
225225
std::size_t gdim = mesh.geometry().dim();
226-
226+
227227
std::vector<std::size_t> columns;
228228
std::vector<double> values;
229-
229+
230230
// Fill the values
231231
for(CellIterator cell(mesh); !cell.end(); ++cell)
232232
{
233233
auto dg_dofs = DG0_dofmap->cell_dofs(cell->index());
234234
// There is only one DG0 dof per cell
235235
dolfin::la_index cell_dof = dg_dofs[0];
236-
236+
237237
Point cell_mp = cell->midpoint();
238238
double cell_volume = cell->volume();
239-
std::size_t local_facet_index = 0;
240-
239+
std::size_t local_facet_index = 0;
240+
241241
auto cr_dofs = CR1_dofmap->cell_dofs(cell->index());
242-
242+
243243
for(FacetIterator facet(*cell); !facet.end(); ++facet)
244244
{
245245
double facet_measure=0;
@@ -248,23 +248,23 @@ void cr_divergence_matrix(GenericMatrix& M, GenericMatrix& A,
248248
else if(tdim == 3)
249249
facet_measure = Face(mesh, facet->index()).area();
250250
// Tdim 1 will not happen because CR is not defined there
251-
251+
252252
Point facet_normal = facet->normal();
253253

254254
// Flip the normal if it is not outer already
255255
Point facet_mp = facet->midpoint();
256256
double sign = (facet_normal.dot(facet_mp - cell_mp) > 0.0) ? 1.0 : -1.0;
257257
facet_normal *= (sign*facet_measure/cell_volume);
258-
258+
259259
// Dofs of CR on the facet, local order
260260
std::vector<std::size_t> facet_dofs;
261261
CR1_dofmap->tabulate_facet_dofs(facet_dofs, local_facet_index);
262-
262+
263263
for (std::size_t j = 0 ; j < facet_dofs.size(); j++)
264-
{
264+
{
265265
columns.push_back(cr_dofs[facet_dofs[j]]);
266266
values.push_back(facet_normal[j]);
267-
}
267+
}
268268
local_facet_index += 1;
269269
}
270270
M.setrow(cell_dof, columns, values);
@@ -275,7 +275,6 @@ void cr_divergence_matrix(GenericMatrix& M, GenericMatrix& A,
275275
//std::shared_ptr<GenericMatrix> Cp = MatMatMult(M, A);
276276
//return Cp;
277277
}
278-
279278
namespace py = pybind11;
280279

281280
PYBIND11_MODULE(cr_divergence, m)
@@ -299,4 +298,3 @@ PYBIND11_MODULE(cr_divergence, m)
299298
cr_divergence_matrix(M, A, _U, _V);
300299
});
301300
}
302-

‎fenicstools/fem/gradient_weight.cpp‎

Lines changed: 41 additions & 42 deletions
Original file line numberDiff line numberDiff line change
@@ -25,27 +25,27 @@ using namespace dolfin;
2525

2626
void compute_weight(Function& DG)
2727
{ // Compute weights for averaging with neighboring cells
28-
28+
2929
// Get the mesh, element and dofmap
30-
std::shared_ptr<const FunctionSpace> V = DG.function_space();
30+
std::shared_ptr<const FunctionSpace> V = DG.function_space();
3131
std::shared_ptr<const Mesh> mesh = V->mesh();
3232
std::shared_ptr<const FiniteElement> element = V->element();
3333
std::shared_ptr<const GenericDofMap> dofmap_u = V->dofmap();
34-
34+
3535
// Allocate storage for weights on one cell
36-
std::vector<double> ws(element->space_dimension());
37-
36+
std::vector<double> ws(element->space_dimension());
37+
3838
// Compute weights
39-
GenericVector& dg_vector = *DG.vector();
39+
GenericVector& dg_vector = *DG.vector();
4040
dg_vector.zero();
4141
for (CellIterator cell(*mesh, "all"); !cell.end(); ++cell)
4242
{
4343
auto dofs = dofmap_u->cell_dofs(cell->index());
44-
44+
4545
std::fill(ws.begin(), ws.end(), 1./cell->volume());
46-
dg_vector.add_local(ws.data(), dofs.size(), dofs.data());
47-
}
48-
dg_vector.apply("insert");
46+
dg_vector.add_local(ws.data(), dofs.size(), dofs.data());
47+
}
48+
dg_vector.apply("insert");
4949
}
5050

5151
std::size_t dof_owner(std::vector<std::vector<std::size_t> > all_ranges,
@@ -66,7 +66,7 @@ void compute_DG0_to_CG_weight_matrix(GenericMatrix& A, Function& DG)
6666
std::vector<double> values;
6767
std::vector<std::vector<std::size_t> > allcolumns;
6868
std::vector<std::vector<double> > allvalues;
69-
69+
7070
const std::pair<std::size_t, std::size_t> row_range = A.local_range(0);
7171
const std::size_t m = row_range.second - row_range.first;
7272
GenericVector& weight = *DG.vector();
@@ -75,26 +75,26 @@ void compute_DG0_to_CG_weight_matrix(GenericMatrix& A, Function& DG)
7575
weight_range_vec[0] = weight_range.first;
7676
weight_range_vec[1] = weight_range.second;
7777
int dm = weight_range.second-weight_range.first;
78-
78+
7979
const MPI_Comm mpi_comm = DG.function_space()->mesh()->mpi_comm();
80-
80+
8181
// Communicate local_ranges of weights
8282
std::vector<std::vector<std::size_t> > all_ranges;
8383
dolfin::MPI::all_gather(mpi_comm, weight_range_vec, all_ranges);
84-
84+
8585
// Number of MPI processes
8686
std::size_t num_processes = dolfin::MPI::size(mpi_comm);
8787

8888
// Some weights live on other processes and need to be communicated
8989
// Create list of off-process weights
90-
std::vector<std::vector<std::size_t> > dofs_needed(num_processes);
90+
std::vector<std::vector<std::size_t> > dofs_needed(num_processes);
9191
for (std::size_t row = 0; row < m; row++)
92-
{
92+
{
9393
// Get global row number
9494
const std::size_t global_row = row + row_range.first;
95-
95+
9696
A.getrow(global_row, columns, values);
97-
97+
9898
for (std::size_t i = 0; i < columns.size(); i++)
9999
{
100100
std::size_t dof = columns[i];
@@ -109,14 +109,14 @@ void compute_DG0_to_CG_weight_matrix(GenericMatrix& A, Function& DG)
109109
// Communicate to all which weights are needed by the process
110110
std::vector<std::vector<std::size_t> > dofs_needed_recv;
111111
dolfin::MPI::all_to_all(mpi_comm, dofs_needed, dofs_needed_recv);
112-
112+
113113
// Fetch the weights that must be communicated
114-
std::vector<std::vector<double> > weights_to_send(num_processes);
114+
std::vector<std::vector<double> > weights_to_send(num_processes);
115115
for (std::size_t p = 0; p < num_processes; p++)
116116
{
117-
if (p == dolfin::MPI::rank(mpi_comm))
117+
if (p == dolfin::MPI::rank(mpi_comm))
118118
continue;
119-
119+
120120
std::vector<std::size_t> dofs = dofs_needed_recv[p];
121121
std::map<std::size_t, double> send_weights;
122122
for (std::size_t k = 0; k < dofs.size(); k++)
@@ -126,25 +126,25 @@ void compute_DG0_to_CG_weight_matrix(GenericMatrix& A, Function& DG)
126126
}
127127
std::vector<std::vector<double> > weights_to_send_recv;
128128
dolfin::MPI::all_to_all(mpi_comm, weights_to_send, weights_to_send_recv);
129-
129+
130130
// Create a map for looking up received weights
131131
std::map<std::size_t, double> received_weights;
132132
for (std::size_t p = 0; p < num_processes; p++)
133133
{
134134
if (p == dolfin::MPI::rank(mpi_comm))
135135
continue;
136-
136+
137137
for (std::size_t k = 0; k < dofs_needed[p].size(); k++)
138138
{
139-
received_weights[dofs_needed[p][k]] = weights_to_send_recv[p][k];
139+
received_weights[dofs_needed[p][k]] = weights_to_send_recv[p][k];
140140
}
141141
}
142-
142+
143143
for (std::size_t row = 0; row < m; row++)
144-
{
144+
{
145145
// Get global row number
146146
const std::size_t global_row = row + row_range.first;
147-
147+
148148
A.getrow(global_row, columns, values);
149149
for (std::size_t i = 0; i < values.size(); i++)
150150
{
@@ -155,14 +155,14 @@ void compute_DG0_to_CG_weight_matrix(GenericMatrix& A, Function& DG)
155155
}
156156
else
157157
{
158-
values[i] = weight[columns[i]-weight_range.first];
158+
values[i] = weight[columns[i]-weight_range.first];
159159
}
160160
// values[i] = 1./values[i];
161161
}
162-
162+
163163
double s = std::accumulate(values.begin(), values.end(), 0.0);
164164
std::transform(values.begin(), values.end(), values.begin(),
165-
std::bind2nd(std::multiplies<double>(), 1./s));
165+
std::bind2nd(std::multiplies<double>(), 1./s));
166166

167167
for (std::size_t i=0; i<values.size(); i++)
168168
{
@@ -174,26 +174,26 @@ void compute_DG0_to_CG_weight_matrix(GenericMatrix& A, Function& DG)
174174
}
175175
else
176176
{
177-
w = weight[dof-weight_range.first];
178-
}
177+
w = weight[dof-weight_range.first];
178+
}
179179
values[i] = values[i]*w;
180180
// values[i] = values[i]*values[i];
181-
182-
}
183-
181+
182+
}
183+
184184
allvalues.push_back(values);
185185
allcolumns.push_back(columns);
186186
}
187187

188188
for (std::size_t row = 0; row < m; row++)
189-
{
189+
{
190190
// Get global row number
191191
const std::size_t global_row = row + row_range.first;
192-
192+
193193
A.setrow(global_row, allcolumns[row], allvalues[row]);
194194
}
195-
A.apply("insert");
196-
}
195+
A.apply("insert");
196+
}
197197

198198
std::shared_ptr<GenericMatrix> MatTransposeMatMult(const GenericMatrix& A, const GenericMatrix& B)
199199
{
@@ -210,7 +210,7 @@ std::shared_ptr<GenericMatrix> MatMatTransposeMult(const GenericMatrix& A, const
210210
const PETScMatrix* Bp = &as_type<const PETScMatrix>(B);
211211
Mat CC;
212212
PetscErrorCode ierr = MatMatTransposeMult(Ap->mat(), Bp->mat(), MAT_INITIAL_MATRIX, PETSC_DEFAULT, &CC);
213-
return PETScMatrix(CC).copy();
213+
return PETScMatrix(CC).copy();
214214
}
215215
//
216216
std::shared_ptr<GenericMatrix> matMatMult(const GenericMatrix& A, const GenericMatrix& B)
@@ -248,4 +248,3 @@ PYBIND11_MODULE(gradient_weight, m)
248248
return MatTransposeMatMult(A, B);
249249
});
250250
}
251-

0 commit comments

Comments
 (0)