Skip to content

Commit

Permalink
Compatible with fenics 2019.1
Browse files Browse the repository at this point in the history
  • Loading branch information
@mikaem
mikaem committed Oct 11, 2019
1 parent 99930d8 commit ad1ee5e
Show file tree
Hide file tree
Showing 4 changed files with 105 additions and 105 deletions.
56 changes: 27 additions & 29 deletions fenicstools/fem/cr_divergence.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ cfg['library_dirs'] = dolfin_pc['library_dirs']

using namespace dolfin;

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

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

// Get topological dimension so that we know what Facet is
const Mesh mesh = *u.function_space()->mesh();
std::size_t tdim = mesh.topology().dim();
std::size_t tdim = mesh.topology().dim();
// Get the info on length of u and gdim for the dot product
std::size_t gdim = mesh.geometry().dim();
std::size_t udim = u.value_dimension(0); // Rank 1!
Expand All @@ -65,10 +65,10 @@ void cr_divergence(Function& divu, const Function& u)
{
Point cell_mp = cell->midpoint();
double cell_volume = cell->volume();

// Dofs of CR on all facets of the cell, global order
auto facets_dofs = CR1_dofmap->cell_dofs(cell->index());

double cell_integral = 0;
std::size_t local_facet_index = 0;
for(FacetIterator facet(*cell); !facet.end(); ++facet)
Expand All @@ -79,7 +79,7 @@ void cr_divergence(Function& divu, const Function& u)
else if(tdim == 3)
facet_measure = Face(mesh, facet->index()).area();
// Tdim 1 will not happen because CR is not defined there

Point facet_normal = facet->normal();

// Flip the normal if it is not outer already
Expand Down Expand Up @@ -115,7 +115,7 @@ void cr_divergence(Function& divu, const Function& u)
void cr_divergence2(Function& divu, const Function& u,
const FunctionSpace& DGscalar, const FunctionSpace& CRvector)
{
// Divu scalar from u vector
// Divu scalar from u vector
std::size_t u_rank = u.value_rank();
std::size_t divu_rank = divu.value_rank();
if((divu_rank == 0) and (u_rank == 1))
Expand All @@ -139,10 +139,10 @@ void cr_divergence2(Function& divu, const Function& u,
std::shared_ptr<const GenericVector> U = u.vector();
std::shared_ptr<GenericVector> DIVU = divu.vector();

// Get dofmaps
// Get dofmaps
std::shared_ptr<const FunctionSpace> DGvector = divu.function_space();
std::shared_ptr<const FunctionSpace> CRtensor = u.function_space();

std::size_t len_divu = divu.value_dimension(0);
// With scalar U can be extracted only once
std::vector<double> U_values;
Expand All @@ -163,7 +163,7 @@ void cr_divergence2(Function& divu, const Function& u,
CRi_rows = CRvector[i]->dofmap()->dofs();
std::size_t m = CRi_rows.size();
U_i->zero();
U_i->set(U_values.data(), m, CRi_rows.data());
U_i->set(U_values.data(), m, CRi_rows.data());
U_i->apply("insert");
}
// For tensor, U_i represents T_ij, (T_i0, T_i1, T_i2)
Expand All @@ -183,7 +183,7 @@ void cr_divergence2(Function& divu, const Function& u,
std::vector<dolfin::la_index>
CRj_rows = CRvector[j]->dofmap()->dofs();

U_i->set(U_values.data(), m, CRj_rows.data());
U_i->set(U_values.data(), m, CRj_rows.data());
U_i->apply("insert");
}
}
Expand All @@ -197,22 +197,22 @@ void cr_divergence2(Function& divu, const Function& u,
std::size_t m = DGi_rows.size();
std::vector<double> DIVU_i_values;
DIVU_i->get_local(DIVU_i_values);
DIVU->set(DIVU_i_values.data(), m, DGi_rows.data());
DIVU->set(DIVU_i_values.data(), m, DGi_rows.data());
DIVU->apply("insert");
}
}
}
// Base case for all divergence computations.
// Base case for all divergence computations.
// Compute divergence of vector field u.
void cr_divergence_matrix(GenericMatrix& M, GenericMatrix& A,
const FunctionSpace& DGscalar,
const FunctionSpace& DGscalar,
const FunctionSpace& CRvector)
{
std::shared_ptr<const GenericDofMap>
CR1_dofmap = CRvector.dofmap(),
DG0_dofmap = DGscalar.dofmap();

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

// Get topological dimension so that we know what Facet is
const Mesh mesh = *DGscalar.mesh();
std::size_t tdim = mesh.topology().dim();
std::size_t tdim = mesh.topology().dim();
std::size_t gdim = mesh.geometry().dim();

std::vector<std::size_t> columns;
std::vector<double> values;

// Fill the values
for(CellIterator cell(mesh); !cell.end(); ++cell)
{
auto dg_dofs = DG0_dofmap->cell_dofs(cell->index());
// There is only one DG0 dof per cell
dolfin::la_index cell_dof = dg_dofs[0];

Point cell_mp = cell->midpoint();
double cell_volume = cell->volume();
std::size_t local_facet_index = 0;
std::size_t local_facet_index = 0;

auto cr_dofs = CR1_dofmap->cell_dofs(cell->index());

for(FacetIterator facet(*cell); !facet.end(); ++facet)
{
double facet_measure=0;
Expand All @@ -248,23 +248,23 @@ void cr_divergence_matrix(GenericMatrix& M, GenericMatrix& A,
else if(tdim == 3)
facet_measure = Face(mesh, facet->index()).area();
// Tdim 1 will not happen because CR is not defined there

Point facet_normal = facet->normal();

// Flip the normal if it is not outer already
Point facet_mp = facet->midpoint();
double sign = (facet_normal.dot(facet_mp - cell_mp) > 0.0) ? 1.0 : -1.0;
facet_normal *= (sign*facet_measure/cell_volume);

// Dofs of CR on the facet, local order
std::vector<std::size_t> facet_dofs;
CR1_dofmap->tabulate_facet_dofs(facet_dofs, local_facet_index);

for (std::size_t j = 0 ; j < facet_dofs.size(); j++)
{
{
columns.push_back(cr_dofs[facet_dofs[j]]);
values.push_back(facet_normal[j]);
}
}
local_facet_index += 1;
}
M.setrow(cell_dof, columns, values);
Expand All @@ -275,7 +275,6 @@ void cr_divergence_matrix(GenericMatrix& M, GenericMatrix& A,
//std::shared_ptr<GenericMatrix> Cp = MatMatMult(M, A);
//return Cp;
}

namespace py = pybind11;

PYBIND11_MODULE(cr_divergence, m)
Expand All @@ -299,4 +298,3 @@ PYBIND11_MODULE(cr_divergence, m)
cr_divergence_matrix(M, A, _U, _V);
});
}

83 changes: 41 additions & 42 deletions fenicstools/fem/gradient_weight.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,27 +25,27 @@ using namespace dolfin;

void compute_weight(Function& DG)
{ // Compute weights for averaging with neighboring cells

// Get the mesh, element and dofmap
std::shared_ptr<const FunctionSpace> V = DG.function_space();
std::shared_ptr<const FunctionSpace> V = DG.function_space();
std::shared_ptr<const Mesh> mesh = V->mesh();
std::shared_ptr<const FiniteElement> element = V->element();
std::shared_ptr<const GenericDofMap> dofmap_u = V->dofmap();

// Allocate storage for weights on one cell
std::vector<double> ws(element->space_dimension());
std::vector<double> ws(element->space_dimension());

// Compute weights
GenericVector& dg_vector = *DG.vector();
GenericVector& dg_vector = *DG.vector();
dg_vector.zero();
for (CellIterator cell(*mesh, "all"); !cell.end(); ++cell)
{
auto dofs = dofmap_u->cell_dofs(cell->index());

std::fill(ws.begin(), ws.end(), 1./cell->volume());
dg_vector.add_local(ws.data(), dofs.size(), dofs.data());
}
dg_vector.apply("insert");
dg_vector.add_local(ws.data(), dofs.size(), dofs.data());
}
dg_vector.apply("insert");
}

std::size_t dof_owner(std::vector<std::vector<std::size_t> > all_ranges,
Expand All @@ -66,7 +66,7 @@ void compute_DG0_to_CG_weight_matrix(GenericMatrix& A, Function& DG)
std::vector<double> values;
std::vector<std::vector<std::size_t> > allcolumns;
std::vector<std::vector<double> > allvalues;

const std::pair<std::size_t, std::size_t> row_range = A.local_range(0);
const std::size_t m = row_range.second - row_range.first;
GenericVector& weight = *DG.vector();
Expand All @@ -75,26 +75,26 @@ void compute_DG0_to_CG_weight_matrix(GenericMatrix& A, Function& DG)
weight_range_vec[0] = weight_range.first;
weight_range_vec[1] = weight_range.second;
int dm = weight_range.second-weight_range.first;

const MPI_Comm mpi_comm = DG.function_space()->mesh()->mpi_comm();

// Communicate local_ranges of weights
std::vector<std::vector<std::size_t> > all_ranges;
dolfin::MPI::all_gather(mpi_comm, weight_range_vec, all_ranges);

// Number of MPI processes
std::size_t num_processes = dolfin::MPI::size(mpi_comm);

// Some weights live on other processes and need to be communicated
// Create list of off-process weights
std::vector<std::vector<std::size_t> > dofs_needed(num_processes);
std::vector<std::vector<std::size_t> > dofs_needed(num_processes);
for (std::size_t row = 0; row < m; row++)
{
{
// Get global row number
const std::size_t global_row = row + row_range.first;

A.getrow(global_row, columns, values);

for (std::size_t i = 0; i < columns.size(); i++)
{
std::size_t dof = columns[i];
Expand All @@ -109,14 +109,14 @@ void compute_DG0_to_CG_weight_matrix(GenericMatrix& A, Function& DG)
// Communicate to all which weights are needed by the process
std::vector<std::vector<std::size_t> > dofs_needed_recv;
dolfin::MPI::all_to_all(mpi_comm, dofs_needed, dofs_needed_recv);

// Fetch the weights that must be communicated
std::vector<std::vector<double> > weights_to_send(num_processes);
std::vector<std::vector<double> > weights_to_send(num_processes);
for (std::size_t p = 0; p < num_processes; p++)
{
if (p == dolfin::MPI::rank(mpi_comm))
if (p == dolfin::MPI::rank(mpi_comm))
continue;

std::vector<std::size_t> dofs = dofs_needed_recv[p];
std::map<std::size_t, double> send_weights;
for (std::size_t k = 0; k < dofs.size(); k++)
Expand All @@ -126,25 +126,25 @@ void compute_DG0_to_CG_weight_matrix(GenericMatrix& A, Function& DG)
}
std::vector<std::vector<double> > weights_to_send_recv;
dolfin::MPI::all_to_all(mpi_comm, weights_to_send, weights_to_send_recv);

// Create a map for looking up received weights
std::map<std::size_t, double> received_weights;
for (std::size_t p = 0; p < num_processes; p++)
{
if (p == dolfin::MPI::rank(mpi_comm))
continue;

for (std::size_t k = 0; k < dofs_needed[p].size(); k++)
{
received_weights[dofs_needed[p][k]] = weights_to_send_recv[p][k];
received_weights[dofs_needed[p][k]] = weights_to_send_recv[p][k];
}
}

for (std::size_t row = 0; row < m; row++)
{
{
// Get global row number
const std::size_t global_row = row + row_range.first;

A.getrow(global_row, columns, values);
for (std::size_t i = 0; i < values.size(); i++)
{
Expand All @@ -155,14 +155,14 @@ void compute_DG0_to_CG_weight_matrix(GenericMatrix& A, Function& DG)
}
else
{
values[i] = weight[columns[i]-weight_range.first];
values[i] = weight[columns[i]-weight_range.first];
}
// values[i] = 1./values[i];
}

double s = std::accumulate(values.begin(), values.end(), 0.0);
std::transform(values.begin(), values.end(), values.begin(),
std::bind2nd(std::multiplies<double>(), 1./s));
std::bind2nd(std::multiplies<double>(), 1./s));

for (std::size_t i=0; i<values.size(); i++)
{
Expand All @@ -174,26 +174,26 @@ void compute_DG0_to_CG_weight_matrix(GenericMatrix& A, Function& DG)
}
else
{
w = weight[dof-weight_range.first];
}
w = weight[dof-weight_range.first];
}
values[i] = values[i]*w;
// values[i] = values[i]*values[i];
}

}

allvalues.push_back(values);
allcolumns.push_back(columns);
}

for (std::size_t row = 0; row < m; row++)
{
{
// Get global row number
const std::size_t global_row = row + row_range.first;

A.setrow(global_row, allcolumns[row], allvalues[row]);
}
A.apply("insert");
}
A.apply("insert");
}

std::shared_ptr<GenericMatrix> MatTransposeMatMult(const GenericMatrix& A, const GenericMatrix& B)
{
Expand All @@ -210,7 +210,7 @@ std::shared_ptr<GenericMatrix> MatMatTransposeMult(const GenericMatrix& A, const
const PETScMatrix* Bp = &as_type<const PETScMatrix>(B);
Mat CC;
PetscErrorCode ierr = MatMatTransposeMult(Ap->mat(), Bp->mat(), MAT_INITIAL_MATRIX, PETSC_DEFAULT, &CC);
return PETScMatrix(CC).copy();
return PETScMatrix(CC).copy();
}
//
std::shared_ptr<GenericMatrix> matMatMult(const GenericMatrix& A, const GenericMatrix& B)
Expand Down Expand Up @@ -248,4 +248,3 @@ PYBIND11_MODULE(gradient_weight, m)
return MatTransposeMatMult(A, B);
});
}

Loading

0 comments on commit ad1ee5e

Please sign in to comment.