diff --git a/demo/Probe.py b/demo/Probe.py index 23b1e99..0d1968f 100644 --- a/demo/Probe.py +++ b/demo/Probe.py @@ -17,11 +17,10 @@ p = Probes(x.flatten(), W) x = x*0.9 p.add_positions(x.flatten(), W) -print 'hei' for i in range(6): p(w0) -print 'hei' print p.array(2, "testarray") # dump snapshot 2 print p.array(filename="testarray") # dump all snapshots print p.dump("testarray") + diff --git a/fenicstools/WeightedGradient.py b/fenicstools/WeightedGradient.py index b4cedf1..ef584f7 100644 --- a/fenicstools/WeightedGradient.py +++ b/fenicstools/WeightedGradient.py @@ -3,7 +3,7 @@ __copyright__ = "Copyright (C) 2013 " + __author__ __license__ = "GNU Lesser GPL version 3 or any later version" import os, inspect -from dolfin import compile_extension_module, Function, FunctionSpace, assemble, TrialFunction, TestFunction, dx, Matrix +from dolfin import info, compile_extension_module, Function, FunctionSpace, assemble, TrialFunction, TestFunction, dx, Matrix fem_folder = os.path.abspath(os.path.join(inspect.getfile(inspect.currentframe()), "../fem")) gradient_code = open(os.path.join(fem_folder, 'gradient_weight.cpp'), 'r').read() @@ -60,6 +60,8 @@ def weighted_gradient_matrix(mesh, i, family='CG', degree=1, constrained_domain= CC.append(Cp) return CC else: - dP = assemble(TrialFunction(S).dx(i)*TestFunction(DG)*dx) + dP = assemble(TrialFunction(S).dx(i)*TestFunction(DG)*dx) Cp = compiled_gradient_module.compute_weighted_gradient_matrix(G, dP, dg) + #info(G, True) + #info(dP, True) return Cp diff --git a/fenicstools/fem/gradient_weight.cpp b/fenicstools/fem/gradient_weight.cpp index f49a23b..e64039e 100644 --- a/fenicstools/fem/gradient_weight.cpp +++ b/fenicstools/fem/gradient_weight.cpp @@ -26,15 +26,15 @@ namespace dolfin // Compute weights GenericVector& dg_vector = *DG.vector(); dg_vector.zero(); - for (CellIterator cell(*mesh); !cell.end(); ++cell) + for (CellIterator cell(*mesh, "all"); !cell.end(); ++cell) { const std::vector& dofs = dofmap_u->cell_dofs(cell->index()); std::fill(ws.begin(), ws.end(), 1./cell->volume()); - dg_vector.add(ws.data(), dofs.size(), dofs.data()); + dg_vector.add_local(ws.data(), dofs.size(), dofs.data()); } - dg_vector.apply("insert"); + dg_vector.apply("insert"); } std::size_t dof_owner(std::vector > all_ranges, @@ -55,7 +55,7 @@ namespace dolfin std::vector values; std::vector > allcolumns; std::vector > allvalues; - + const std::pair row_range = A.local_range(0); const std::size_t m = row_range.second - row_range.first; GenericVector& weight = *DG.vector(); @@ -66,7 +66,7 @@ namespace dolfin 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 > all_ranges; MPI::all_gather(mpi_comm, weight_range_vec, all_ranges); @@ -94,7 +94,7 @@ namespace dolfin } } } - + // Communicate to all which weights are needed by the process std::vector > dofs_needed_recv; MPI::all_to_all(mpi_comm, dofs_needed, dofs_needed_recv); @@ -110,7 +110,7 @@ namespace dolfin std::map send_weights; for (std::size_t k = 0; k < dofs.size(); k++) { - weights_to_send[p].push_back(weight[dofs[k]]); + weights_to_send[p].push_back(weight[dofs[k]-weight_range.first]); } } std::vector > weights_to_send_recv; @@ -124,7 +124,9 @@ namespace dolfin 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++) @@ -142,7 +144,7 @@ namespace dolfin } else { - values[i] = weight[columns[i]]; + values[i] = weight[columns[i]-weight_range.first]; } // values[i] = 1./values[i]; } @@ -161,7 +163,7 @@ namespace dolfin } else { - w = weight[dof]; + w = weight[dof-weight_range.first]; } values[i] = values[i]*w; // values[i] = values[i]*values[i]; @@ -171,6 +173,7 @@ namespace dolfin allvalues.push_back(values); allcolumns.push_back(columns); } + for (std::size_t row = 0; row < m; row++) { // Get global row number @@ -178,7 +181,7 @@ namespace dolfin A.setrow(global_row, allcolumns[row], allvalues[row]); } - A.apply("insert"); + A.apply("insert"); } std::shared_ptr MatTransposeMatMult(GenericMatrix& A, GenericMatrix& B) @@ -208,6 +211,7 @@ namespace dolfin Mat CC; PetscErrorCode ierr = MatMatMult(Ap->mat(), Bp->mat(), MAT_INITIAL_MATRIX, PETSC_DEFAULT, &CC); dolfin::PETScMatrix CCC = PETScMatrix(CC); + CCC.apply("insert"); return CCC.copy(); } diff --git a/tests/test_StatisticsProbes.py b/tests/test_StatisticsProbes.py index b2387fd..eabc38b 100644 --- a/tests/test_StatisticsProbes.py +++ b/tests/test_StatisticsProbes.py @@ -1,7 +1,7 @@ import nose from dolfin import FunctionSpace, UnitCubeMesh, UnitSquareMesh, interpolate, \ - Expression, VectorFunctionSpace + Expression, VectorFunctionSpace, MPI, mpi_comm_world from fenicstools import * from numpy import array @@ -17,22 +17,10 @@ def test_StatisticsProbes_segregated_2D(): for i in range(5): probes(u0, v0) - id0, p = probes[0] - if len(p) > 0: - assert p.number_of_evaluations() == 5 - assert p.value_size() == 5 - - mean = p.mean() - var = p.variance() - nose.tools.assert_almost_equal(p[0][0], 2.5) - nose.tools.assert_almost_equal(p[0][4], 0.625) - nose.tools.assert_almost_equal(p[1][0], 0.5) - nose.tools.assert_almost_equal(p[1][1], 0.25) - nose.tools.assert_almost_equal(mean[0], 0.5) - nose.tools.assert_almost_equal(mean[1], 0.25) - nose.tools.assert_almost_equal(var[0], 0.25) - nose.tools.assert_almost_equal(var[1], 0.0625) - nose.tools.assert_almost_equal(var[2], 0.125) + p = probes.array() + if MPI.rank(mpi_comm_world()) == 0: + nose.tools.assert_almost_equal(p[0,0], 2.5) + nose.tools.assert_almost_equal(p[0,4], 0.625) def test_StatisticsProbes_segregated_3D(): mesh = UnitCubeMesh(4, 4, 4) @@ -47,26 +35,10 @@ def test_StatisticsProbes_segregated_3D(): for i in range(5): probes(u0, v0, w0) - id0, p = probes[0] - if len(p) > 0: - assert p.number_of_evaluations() == 5 - assert p.value_size() == 9 - - mean = p.mean() - var = p.variance() - nose.tools.assert_almost_equal(p[0][0], 2.5) - nose.tools.assert_almost_equal(p[0][4], 0.3125) - nose.tools.assert_almost_equal(p[1][0], 0.5) - nose.tools.assert_almost_equal(p[1][1], 0.25) - nose.tools.assert_almost_equal(p[1][2], 0.25) - nose.tools.assert_almost_equal(mean[0], 0.5) - nose.tools.assert_almost_equal(mean[1], 0.25) - nose.tools.assert_almost_equal(mean[2], 0.25) - nose.tools.assert_almost_equal(var[0], 0.25) - nose.tools.assert_almost_equal(var[1], 0.0625) - nose.tools.assert_almost_equal(var[2], 0.0625) - nose.tools.assert_almost_equal(var[3], 0.125) - nose.tools.assert_almost_equal(var[4], 0.125) + p = probes.array() + if MPI.rank(mpi_comm_world()) == 0: + nose.tools.assert_almost_equal(p[0,0], 2.5) + nose.tools.assert_almost_equal(p[0,4], 0.3125) def test_StatisticsProbes_vector_2D(): mesh = UnitSquareMesh(4, 4) @@ -79,22 +51,10 @@ def test_StatisticsProbes_vector_2D(): for i in range(5): probes(u0) - id0, p = probes[0] - if len(p) > 0: - assert p.number_of_evaluations() == 5 - assert p.value_size() == 5 - - mean = p.mean() - var = p.variance() - nose.tools.assert_almost_equal(p[0][0], 2.5) - nose.tools.assert_almost_equal(p[0][4], 0.625) - nose.tools.assert_almost_equal(p[1][0], 0.5) - nose.tools.assert_almost_equal(p[1][1], 0.25) - nose.tools.assert_almost_equal(mean[0], 0.5) - nose.tools.assert_almost_equal(mean[1], 0.25) - nose.tools.assert_almost_equal(var[0], 0.25) - nose.tools.assert_almost_equal(var[1], 0.0625) - nose.tools.assert_almost_equal(var[2], 0.125) + p = probes.array() + if MPI.rank(mpi_comm_world()) == 0: + nose.tools.assert_almost_equal(p[0,0], 2.5) + nose.tools.assert_almost_equal(p[0,4], 0.625) def test_StatisticsProbes_vector_3D(): mesh = UnitCubeMesh(4, 4, 4) @@ -107,26 +67,11 @@ def test_StatisticsProbes_vector_3D(): for i in range(5): probes(u0) - id0, p = probes[0] - if len(p) > 0: - assert p.number_of_evaluations() == 5 - assert p.value_size() == 9 + p = probes.array() + if MPI.rank(mpi_comm_world()) == 0: - mean = p.mean() - var = p.variance() - nose.tools.assert_almost_equal(p[0][0], 2.5) - nose.tools.assert_almost_equal(p[0][4], 0.3125) - nose.tools.assert_almost_equal(p[1][0], 0.5) - nose.tools.assert_almost_equal(p[1][1], 0.25) - nose.tools.assert_almost_equal(p[1][2], 0.25) - nose.tools.assert_almost_equal(mean[0], 0.5) - nose.tools.assert_almost_equal(mean[1], 0.25) - nose.tools.assert_almost_equal(mean[2], 0.25) - nose.tools.assert_almost_equal(var[0], 0.25) - nose.tools.assert_almost_equal(var[1], 0.0625) - nose.tools.assert_almost_equal(var[2], 0.0625) - nose.tools.assert_almost_equal(var[3], 0.125) - nose.tools.assert_almost_equal(var[4], 0.125) + nose.tools.assert_almost_equal(p[0,0], 2.5) + nose.tools.assert_almost_equal(p[0,4], 0.3125) -if __name__ == '__main__': - nose.run(defaultTest=__name__) \ No newline at end of file +#if __name__ == '__main__': + #nose.run(defaultTest=__name__) \ No newline at end of file diff --git a/tests/test_WeightedGradient.py b/tests/test_WeightedGradient.py index 76a229b..4f48f2d 100644 --- a/tests/test_WeightedGradient.py +++ b/tests/test_WeightedGradient.py @@ -30,5 +30,5 @@ def test_WeightedGradient(): du.vector()[:] = wx[2] * u.vector() nose.tools.assert_almost_equal(du.vector().min(), 3.) -#if __name__ == '__main__': - #nose.run(defaultTest=__name__) \ No newline at end of file +if __name__ == '__main__': + nose.run(defaultTest=__name__) \ No newline at end of file