Merge pull request #3443 from grmnptr/piecewise-divide

Add pointwise_divide to PetscVector

Author: lindsayad

include/numerics/distributed_vector.h

@@ -226,6 +226,9 @@ class DistributedVector final : public NumericVector<T>
   virtual void pointwise_mult (const NumericVector<T> & vec1,
                                const NumericVector<T> & vec2) override;
 
+  virtual void pointwise_divide (const NumericVector<T> & vec1,
+                                 const NumericVector<T> & vec2) override;
+
   virtual void swap (NumericVector<T> & v) override;
 
   virtual std::size_t max_allowed_id() const override;

include/numerics/eigen_sparse_vector.h

@@ -228,6 +228,9 @@ class EigenSparseVector final : public NumericVector<T>
   virtual void pointwise_mult (const NumericVector<T> & vec1,
                                const NumericVector<T> & vec2) override;
 
+  virtual void pointwise_divide (const NumericVector<T> & vec1,
+                                 const NumericVector<T> & vec2) override;
+
   virtual void swap (NumericVector<T> & v) override;
 
   virtual std::size_t max_allowed_id() const override;

include/numerics/laspack_vector.h

@@ -226,6 +226,9 @@ class LaspackVector final : public NumericVector<T>
   virtual void pointwise_mult (const NumericVector<T> & vec1,
                                const NumericVector<T> & vec2) override;
 
+  virtual void pointwise_divide (const NumericVector<T> & vec1,
+                                 const NumericVector<T> & vec2) override;
+
   virtual void swap (NumericVector<T> & v) override;
 
   virtual std::size_t max_allowed_id() const override;

include/numerics/numeric_vector.h

@@ -671,6 +671,14 @@ class NumericVector : public ReferenceCountedObject<NumericVector<T>>,
   virtual void pointwise_mult (const NumericVector<T> & vec1,
                                const NumericVector<T> & vec2) = 0;
 
+  /**
+   * Computes \f$ u_i \leftarrow \frac{v_{1,i}}{v_{2,i}} \f$ (summation not implied)
+   * i.e. the pointwise (component-wise) division of \p vec1 and
+   * \p vec2, and stores the result in \p *this.
+   */
+  virtual void pointwise_divide (const NumericVector<T> & vec1,
+                                 const NumericVector<T> & vec2) = 0;
+
   /**
    * Prints the local contents of the vector, by default to
    * libMesh::out

include/numerics/petsc_vector.h

@@ -321,6 +321,9 @@ class PetscVector final : public NumericVector<T>
   virtual void pointwise_mult (const NumericVector<T> & vec1,
                                const NumericVector<T> & vec2) override;
 
+  virtual void pointwise_divide (const NumericVector<T> & vec1,
+                                 const NumericVector<T> & vec2) override;
+
   virtual void print_matlab(const std::string & name = "") const override;
 
   virtual void create_subvector(NumericVector<T> & subvector,

src/numerics/distributed_vector.C

@@ -619,7 +619,12 @@ void DistributedVector<T>::pointwise_mult (const NumericVector<T> &,
   libmesh_not_implemented();
 }
 
-
+template <typename T>
+void DistributedVector<T>::pointwise_divide (const NumericVector<T> &,
+                                             const NumericVector<T> &)
+{
+  libmesh_not_implemented();
+}
 
 //--------------------------------------------------------------
 // Explicit instantiations

src/numerics/eigen_sparse_vector.C

@@ -418,6 +418,12 @@ void EigenSparseVector<T>::pointwise_mult (const NumericVector<T> & /*vec1*/,
   libmesh_not_implemented();
 }
 
+template <typename T>
+void EigenSparseVector<T>::pointwise_divide (const NumericVector<T> & /*vec1*/,
+                                             const NumericVector<T> & /*vec2*/)
+{
+  libmesh_not_implemented();
+}
 
 
 template <typename T>

src/numerics/laspack_vector.C

@@ -482,7 +482,12 @@ void LaspackVector<T>::pointwise_mult (const NumericVector<T> & /*vec1*/,
   libmesh_not_implemented();
 }
 
-
+template <typename T>
+void LaspackVector<T>::pointwise_divide (const NumericVector<T> & /*vec1*/,
+                                         const NumericVector<T> & /*vec2*/)
+{
+  libmesh_not_implemented();
+}
 
 template <typename T>
 Real LaspackVector<T>::max() const

src/numerics/petsc_vector.C

@@ -1022,7 +1022,24 @@ void PetscVector<T>::pointwise_mult (const NumericVector<T> & vec1,
   this->_is_closed = true;
 }
 
+template <typename T>
+void PetscVector<T>::pointwise_divide (const NumericVector<T> & vec1,
+                                       const NumericVector<T> & vec2)
+{
+  this->_restore_array();
+
+  // Convert arguments to PetscVector*.
+  const PetscVector<T> * const vec1_petsc = cast_ptr<const PetscVector<T> *>(&vec1);
+  const PetscVector<T> * const vec2_petsc = cast_ptr<const PetscVector<T> *>(&vec2);
+
+  // Call PETSc function.
+  PetscErrorCode ierr = VecPointwiseDivide(_vec, vec1_petsc->vec(), vec2_petsc->vec());
+  LIBMESH_CHKERR(ierr);
+  if (this->type() == GHOSTED)
+    VecGhostUpdateBeginEnd(this->comm(), _vec, INSERT_VALUES, SCATTER_FORWARD);
 
+  this->_is_closed = true;
+}
 
 template <typename T>
 void PetscVector<T>::print_matlab (const std::string & name) const

tests/numerics/petsc_vector_test.C

@@ -23,6 +23,8 @@ public:
 
   CPPUNIT_TEST( testGetArray );
 
+  CPPUNIT_TEST( testPetscOperations );
+
   CPPUNIT_TEST_SUITE_END();
 
   void testGetArray()
@@ -77,6 +79,41 @@ public:
     v.restore_array();
   }
 
+  void testPetscOperations()
+  {
+    PetscVector<Number> v1(*my_comm, global_size, local_size);
+    auto v2 = v1.clone();
+
+    const libMesh::dof_id_type
+      first = v1.first_local_index(),
+      last  = v1.last_local_index();
+
+    for (libMesh::dof_id_type n=first; n != last; n++)
+    {
+      v1.set (n, static_cast<libMesh::Number>(n+1));
+      v2->set (n, static_cast<libMesh::Number>(2*(n+1)));
+    }
+    v1.close();
+    v2->close();
+
+    auto v_working_ptr = v1.clone();
+    auto & v_working = *v_working_ptr;
+
+    v_working.pointwise_mult(v1, *v2);
+
+    for (libMesh::dof_id_type n=first; n != last; n++)
+      LIBMESH_ASSERT_FP_EQUAL(libMesh::libmesh_real(v_working(n)),
+                              libMesh::Real((n+1)*2*(n+1)),
+                              libMesh::TOLERANCE*libMesh::TOLERANCE);
+
+    v_working.pointwise_divide(v1, *v2);
+
+    for (libMesh::dof_id_type n=first; n != last; n++)
+      LIBMESH_ASSERT_FP_EQUAL(libMesh::libmesh_real(v_working(n)),
+                              libMesh::Real(0.5),
+                              libMesh::TOLERANCE*libMesh::TOLERANCE);
+  }
+
 };
 
 CPPUNIT_TEST_SUITE_REGISTRATION( PetscVectorTest );