diff --git a/src/DCPSE/DCPSE_op/DCPSE_surface_op.hpp b/src/DCPSE/DCPSE_op/DCPSE_surface_op.hpp
index 030c78730c006aa59fa6d2dd8dd27994cbe31e9b..eed13a17b1e8be32ffad5351f7886efaff1f2490 100644
--- a/src/DCPSE/DCPSE_op/DCPSE_surface_op.hpp
+++ b/src/DCPSE/DCPSE_op/DCPSE_surface_op.hpp
@@ -62,6 +62,27 @@ public:
}
+ /*! \brief Method for Saving the DCPSE Operator.
+ *
+ * \param parts particle set
+ * \param file name for data to be saved.
+ */
+ template<typename particles_type>
+ void save(particles_type &particles, const std::string &file) {
+ auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
+ dcpse_temp->save(file);
+ }
+ /*! \brief Method for Loading the DCPSE Operator.
+ *
+ * \param parts particle set
+ * \param file name for data to be loaded from.
+ */
+ template<typename particles_type>
+ void load(particles_type &particles, const std::string &file) {
+ auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
+ dcpse_temp->load(file);
+ }
+
/*! \brief Method for Updating the DCPSE Operator by recomputing DCPSE Kernels.
*
*
@@ -131,6 +152,27 @@ public:
}
+
+ /*! \brief Method for Saving the DCPSE Operator.
+ *
+ * \param parts particle set
+ * \param file name for data to be saved.
+ */
+ template<typename particles_type>
+ void save(particles_type &particles, const std::string &file) {
+ auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
+ dcpse_temp->save(file);
+ }
+ /*! \brief Method for Loading the DCPSE Operator.
+ *
+ * \param parts particle set
+ * \param file name for data to be loaded from.
+ */
+ template<typename particles_type>
+ void load(particles_type &particles, const std::string &file) {
+ auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
+ dcpse_temp->load(file);
+ }
/*! \brief Method for Updating the DCPSE Operator by recomputing DCPSE Kernels.
*
*
@@ -201,6 +243,26 @@ public:
}
+ /*! \brief Method for Saving the DCPSE Operator.
+ *
+ * \param parts particle set
+ * \param file name for data to be saved.
+ */
+ template<typename particles_type>
+ void save(particles_type &particles, const std::string &file) {
+ auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
+ dcpse_temp->save(file);
+ }
+ /*! \brief Method for Loading the DCPSE Operator.
+ *
+ * \param parts particle set
+ * \param file name for data to be loaded from.
+ */
+ template<typename particles_type>
+ void load(particles_type &particles, const std::string &file) {
+ auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
+ dcpse_temp->load(file);
+ }
/*! \brief Method for Updating the DCPSE Operator by recomputing DCPSE Kernels.
*
*
@@ -273,6 +335,26 @@ public:
}
+ /*! \brief Method for Saving the DCPSE Operator.
+ *
+ * \param parts particle set
+ * \param file name for data to be saved.
+ */
+ template<typename particles_type>
+ void save(particles_type &particles, const std::string &file) {
+ auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
+ dcpse_temp->save(file);
+ }
+ /*! \brief Method for Loading the DCPSE Operator.
+ *
+ * \param parts particle set
+ * \param file name for data to be loaded from.
+ */
+ template<typename particles_type>
+ void load(particles_type &particles, const std::string &file) {
+ auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
+ dcpse_temp->load(file);
+ }
/*! \brief Method for Updating the DCPSE Operator by recomputing DCPSE Kernels.
*
*
@@ -343,6 +425,26 @@ public:
}
+ /*! \brief Method for Saving the DCPSE Operator.
+ *
+ * \param parts particle set
+ * \param file name for data to be saved.
+ */
+ template<typename particles_type>
+ void save(particles_type &particles, const std::string &file) {
+ auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
+ dcpse_temp->save(file);
+ }
+ /*! \brief Method for Loading the DCPSE Operator.
+ *
+ * \param parts particle set
+ * \param file name for data to be loaded from.
+ */
+ template<typename particles_type>
+ void load(particles_type &particles, const std::string &file) {
+ auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
+ dcpse_temp->load(file);
+ }
/*! \brief Method for Updating the DCPSE Operator by recomputing DCPSE Kernels.
*
*
@@ -414,6 +516,26 @@ public:
}
+ /*! \brief Method for Saving the DCPSE Operator.
+ *
+ * \param parts particle set
+ * \param file name for data to be saved.
+ */
+ template<typename particles_type>
+ void save(particles_type &particles, const std::string &file) {
+ auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
+ dcpse_temp->save(file);
+ }
+ /*! \brief Method for Loading the DCPSE Operator.
+ *
+ * \param parts particle set
+ * \param file name for data to be loaded from.
+ */
+ template<typename particles_type>
+ void load(particles_type &particles, const std::string &file) {
+ auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
+ dcpse_temp->load(file);
+ }
/*! \brief Method for Updating the DCPSE Operator by recomputing DCPSE Kernels.
*
*
@@ -484,6 +606,26 @@ public:
}
+ /*! \brief Method for Saving the DCPSE Operator.
+ *
+ * \param parts particle set
+ * \param file name for data to be saved.
+ */
+ template<typename particles_type>
+ void save(particles_type &particles, const std::string &file) {
+ auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
+ dcpse_temp->save(file);
+ }
+ /*! \brief Method for Loading the DCPSE Operator.
+ *
+ * \param parts particle set
+ * \param file name for data to be loaded from.
+ */
+ template<typename particles_type>
+ void load(particles_type &particles, const std::string &file) {
+ auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
+ dcpse_temp->load(file);
+ }
/*! \brief Method for Updating the DCPSE Operator by recomputing DCPSE Kernels.
*
*
@@ -556,6 +698,26 @@ public:
}
+ /*! \brief Method for Saving the DCPSE Operator.
+ *
+ * \param parts particle set
+ * \param file name for data to be saved.
+ */
+ template<typename particles_type>
+ void save(particles_type &particles, const std::string &file) {
+ auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
+ dcpse_temp->save(file);
+ }
+ /*! \brief Method for Loading the DCPSE Operator.
+ *
+ * \param parts particle set
+ * \param file name for data to be loaded from.
+ */
+ template<typename particles_type>
+ void load(particles_type &particles, const std::string &file) {
+ auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
+ dcpse_temp->load(file);
+ }
/*! \brief Method for Updating the DCPSE Operator by recomputing DCPSE Kernels.
*
*
@@ -628,6 +790,26 @@ public:
}
+ /*! \brief Method for Saving the DCPSE Operator.
+ *
+ * \param parts particle set
+ * \param file name for data to be saved.
+ */
+ template<typename particles_type>
+ void save(particles_type &particles, const std::string &file) {
+ auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
+ dcpse_temp->save(file);
+ }
+ /*! \brief Method for Loading the DCPSE Operator.
+ *
+ * \param parts particle set
+ * \param file name for data to be loaded from.
+ */
+ template<typename particles_type>
+ void load(particles_type &particles, const std::string &file) {
+ auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
+ dcpse_temp->load(file);
+ }
/*! \brief Method for Updating the DCPSE Operator by recomputing DCPSE Kernels.
*
*
@@ -700,6 +882,113 @@ public:
}
+ /*! \brief Method for Saving the DCPSE Operator.
+ *
+ * \param parts particle set
+ * \param file name for data to be saved.
+ */
+ template<typename particles_type>
+ void save(particles_type &particles, const std::string &file) {
+ auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
+ dcpse_temp->save(file);
+ }
+ /*! \brief Method for Loading the DCPSE Operator.
+ *
+ * \param parts particle set
+ * \param file name for data to be loaded from.
+ */
+ template<typename particles_type>
+ void load(particles_type &particles, const std::string &file) {
+ auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
+ dcpse_temp->load(file);
+ }
+ /*! \brief Method for Updating the DCPSE Operator by recomputing DCPSE Kernels.
+ *
+ *
+ * \param parts particle set
+ */
+ template<typename particles_type>
+ void update(particles_type &particles) {
+ auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
+ dcpse_temp->createNormalParticles<NORMAL_ID>(particles);
+ dcpse_temp->initializeUpdate(particles);
+ dcpse_temp->accumulateAndDeleteNormalParticles(particles);
+
+ }
+};
+
+template<unsigned int NORMAL_ID>
+class SurfaceDerivative_G {
+
+ void *dcpse;
+
+public:
+ /*! \brief Class for Creating the DCPSE Operator Dxx and objects and computs DCPSE Kernels.
+ *
+ *
+ * \param parts particle set
+ * \param ord order of convergence of the operator
+ * \param rCut Argument for cell list construction
+ * \param oversampling_factor multiplier to the minimum no. of particles required by the operator in support
+ * \param support_options default:N_particles, Radius can be used to select all particles inside rCut. Overrides oversampling.
+ *
+ * \return Operator Dxx which is a function on Vector_dist_Expressions
+ *
+ */
+ template<typename particles_type>
+ SurfaceDerivative_G(particles_type &parts, unsigned int ord, typename particles_type::stype rCut,typename particles_type::stype nSpacing,
+ const Point<particles_type::dims, unsigned int> &p,support_options opt = support_options::RADIUS) {
+
+ dcpse = new Dcpse<particles_type::dims, particles_type>(parts, p, ord, rCut,nSpacing,value_t<NORMAL_ID>(), opt);
+ }
+
+ template<typename particles_type>
+ void deallocate(particles_type &parts) {
+ delete (Dcpse<particles_type::dims, particles_type> *) dcpse;
+ }
+
+ template<typename operand_type>
+ vector_dist_expression_op<operand_type, Dcpse<operand_type::vtype::dims, typename operand_type::vtype>, VECT_DCPSE>
+ operator()(operand_type arg) {
+ typedef Dcpse<operand_type::vtype::dims, typename operand_type::vtype> dcpse_type;
+
+ return vector_dist_expression_op<operand_type, dcpse_type, VECT_DCPSE>(arg, *(dcpse_type *) dcpse);
+ }
+
+ template<typename particles_type>
+ void checkMomenta(particles_type &particles) {
+ auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
+ dcpse_temp->checkMomenta(particles);
+
+ }
+
+ template<unsigned int prp, typename particles_type>
+ void DrawKernel(particles_type &particles, int k) {
+ auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
+ dcpse_temp->template DrawKernel<prp>(particles, k);
+
+ }
+
+ /*! \brief Method for Saving the DCPSE Operator.
+ *
+ * \param parts particle set
+ * \param file name for data to be saved.
+ */
+ template<typename particles_type>
+ void save(particles_type &particles, const std::string &file) {
+ auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
+ dcpse_temp->save(file);
+ }
+ /*! \brief Method for Loading the DCPSE Operator.
+ *
+ * \param parts particle set
+ * \param file name for data to be loaded from.
+ */
+ template<typename particles_type>
+ void load(particles_type &particles, const std::string &file) {
+ auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
+ dcpse_temp->load(file);
+ }
/*! \brief Method for Updating the DCPSE Operator by recomputing DCPSE Kernels.
*
*
diff --git a/src/DCPSE/DCPSE_op/tests/DCPSE_op_Surface_tests.cpp b/src/DCPSE/DCPSE_op/tests/DCPSE_op_Surface_tests.cpp
index ed225ee05dc2345f1e918bf2067e3944520e8052..3664e8627074a4b2c4380c00f39c0a5480aa32c8 100644
--- a/src/DCPSE/DCPSE_op/tests/DCPSE_op_Surface_tests.cpp
+++ b/src/DCPSE/DCPSE_op/tests/DCPSE_op_Surface_tests.cpp
@@ -645,6 +645,250 @@ BOOST_AUTO_TEST_CASE(dcpse_surface_sphere_old) {
vector_dist<3, double, aggregate<double,double,double,double,double,double,double[3]>> domain(0, box, bc, ghost);
+ //Init_DCPSE(domain)
+ BOOST_TEST_MESSAGE("Init domain...");
+
+ auto it = domain.getGridIterator(sz);
+ while (it.isNext()) {
+ domain.add();
+
+ auto key = it.get();
+ double x = key.get(0) * it.getSpacing(0);
+ domain.getLastPos()[0] = x;
+ double y = key.get(1) * it.getSpacing(1);
+ domain.getLastPos()[1] = y;
+
+ domain.getLastPos()[2] = 0;
+
+ ++it;
+ }
+
+ // Add multi res patch 1
+
+ {
+ const size_t sz2[3] = {40,40,1};
+ Box<3,double> bx({0.25 + it.getSpacing(0)/4.0,0.25 + it.getSpacing(0)/4.0,-0.5},{sz2[0]*it.getSpacing(0)/2.0 + 0.25 + it.getSpacing(0)/4.0, sz2[1]*it.getSpacing(0)/2.0 + 0.25 + it.getSpacing(0)/4.0,0.5});
+ openfpm::vector<size_t> rem;
+
+ auto it = domain.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto k = it.get();
+
+ Point<3,double> xp = domain.getPos(k);
+
+ if (bx.isInside(xp) == true)
+ {
+ rem.add(k.getKey());
+ }
+
+ ++it;
+ }
+
+ domain.remove(rem);
+
+ auto it2 = domain.getGridIterator(sz2);
+ while (it2.isNext()) {
+ domain.add();
+
+ auto key = it2.get();
+ double x = key.get(0) * spacing/2.0 + 0.25 + spacing/4.0;
+ domain.getLastPos()[0] = x;
+ double y = key.get(1) * spacing/2.0 + 0.25 + spacing/4.0;
+ domain.getLastPos()[1] = y;
+ domain.getLastPos()[2] = 0;
+
+ ++it2;
+ }
+ }
+
+ // Add multi res patch 2
+
+ {
+ const size_t sz2[3] = {40,40,1};
+ Box<3,double> bx({0.25 + 21.0*spacing/8.0,0.25 + 21.0*spacing/8.0,-5},{sz2[0]*spacing/4.0 + 0.25 + 21.0*spacing/8.0, sz2[1]*spacing/4.0 + 0.25 + 21*spacing/8.0,5});
+ openfpm::vector<size_t> rem;
+
+ auto it = domain.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto k = it.get();
+
+ Point<3,double> xp = domain.getPos(k);
+
+ if (bx.isInside(xp) == true)
+ {
+ rem.add(k.getKey());
+ }
+
+ ++it;
+ }
+
+ domain.remove(rem);
+
+ auto it2 = domain.getGridIterator(sz2);
+ while (it2.isNext()) {
+ domain.add();
+ auto key = it2.get();
+ double x = key.get(0) * spacing/4.0 + 0.25 + 21*spacing/8.0;
+ domain.getLastPos()[0] = x;
+ double y = key.get(1) * spacing/4.0 + 0.25 + 21*spacing/8.0;
+ domain.getLastPos()[1] = y;
+ domain.getLastPos()[2] = 0;
+
+ ++it2;
+ }
+ }
+
+ ///////////////////////
+
+ BOOST_TEST_MESSAGE("Sync domain across processors...");
+
+ domain.map();
+ domain.ghost_get<0>();
+
+ openfpm::vector<aggregate<int>> bulk;
+ openfpm::vector<aggregate<int>> up_p;
+ openfpm::vector<aggregate<int>> dw_p;
+ openfpm::vector<aggregate<int>> l_p;
+ openfpm::vector<aggregate<int>> r_p;
+ openfpm::vector<aggregate<int>> ref_p;
+
+ auto v = getV<0>(domain);
+ auto RHS=getV<1>(domain);
+ auto sol = getV<2>(domain);
+ auto anasol = getV<3>(domain);
+ auto err = getV<4>(domain);
+ auto DCPSE_sol=getV<5>(domain);
+
+ // Here fill me
+
+ Box<3, double> up({box.getLow(0) - spacing / 2.0, box.getHigh(1) - spacing / 2.0,-5},
+ {box.getHigh(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0,5});
+
+ Box<3, double> down({box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0,-5},
+ {box.getHigh(0) + spacing / 2.0, box.getLow(1) + spacing / 2.0,5});
+
+ Box<3, double> left({box.getLow(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0,-5},
+ {box.getLow(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0,5});
+
+ Box<3, double> right({box.getHigh(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0,-5},
+ {box.getHigh(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0,5});
+
+ openfpm::vector<Box<3, double>> boxes;
+ boxes.add(up);
+ boxes.add(down);
+ boxes.add(left);
+ boxes.add(right);
+
+ // Create a writer and write
+ VTKWriter<openfpm::vector<Box<3, double>>, VECTOR_BOX> vtk_box;
+ vtk_box.add(boxes);
+ //vtk_box.write("vtk_box.vtk");
+
+ auto it2 = domain.getDomainIterator();
+
+ while (it2.isNext()) {
+ auto p = it2.get();
+ Point<3, double> xp = domain.getPos(p);
+ if (up.isInside(xp) == true) {
+ up_p.add();
+ up_p.last().get<0>() = p.getKey();
+ domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+ domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+ } else if (down.isInside(xp) == true) {
+ dw_p.add();
+ dw_p.last().get<0>() = p.getKey();
+ domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+ domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+
+ } else if (left.isInside(xp) == true) {
+ l_p.add();
+ l_p.last().get<0>() = p.getKey();
+ domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+ domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+
+ } else if (right.isInside(xp) == true) {
+ r_p.add();
+ r_p.last().get<0>() = p.getKey();
+ domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+ domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+
+ } else {
+ bulk.add();
+ bulk.last().get<0>() = p.getKey();
+ domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+ domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+ }
+ domain.getProp<6>(p)[0] = 0;
+ domain.getProp<6>(p)[1] = 0;
+ domain.getProp<6>(p)[2] = 1;
+
+ ++it2;
+ }
+
+ domain.ghost_get<1,2,3>();
+ SurfaceDerivative_xx<6> Dxx(domain, 2, rCut,1.0,support_options::ADAPTIVE_SURFACE);
+
+/* v=0;
+ auto itNNN=domain.getDomainIterator();
+ while(itNNN.isNext()){
+ auto p=itNNN.get().getKey();
+ Dxx.DrawKernel<0,decltype(domain)>(domain,p);
+ domain.write_frame("Kernel",p);
+ v=0;
+ ++itNNN;
+ }
+
+*/
+ //Dxx.DrawKernel<5,decltype(domain)>(domain,6161);
+ //domain.write_frame("Kernel",6161);
+
+ SurfaceDerivative_yy<6> Dyy(domain, 2, rCut,1.0,support_options::ADAPTIVE_SURFACE);
+ SurfaceDerivative_zz<6> Dzz(domain, 2, rCut,1.0,support_options::ADAPTIVE_SURFACE);
+
+ Dxx.save(domain,"Sdxx_test");
+ Dyy.save(domain,"Sdyy_test");
+ Dzz.save(domain,"Sdzz_test");
+
+ domain.ghost_get<2>();
+ sol=Dxx(anasol)+Dyy(anasol)+Dzz(anasol);
+ domain.ghost_get<5>();
+
+ double worst1 = 0.0;
+
+ for(int j=0;j<bulk.size();j++)
+ { auto p=bulk.get<0>(j);
+ if (fabs(domain.getProp<1>(p) - domain.getProp<2>(p)) >= worst1) {
+ worst1 = fabs(domain.getProp<1>(p) - domain.getProp<2>(p));
+ }
+ domain.getProp<4>(p) = fabs(domain.getProp<1>(p) - domain.getProp<2>(p));
+
+ }
+ //std::cout << "Maximum Analytic Error: " << worst1 << std::endl;
+ //domain.ghost_get<4>();
+ //domain.write("Robin_anasol");
+ BOOST_REQUIRE(worst1 < 0.03);
+
+ }
+
+
+ BOOST_AUTO_TEST_CASE(dcpse_surface_adaptive_load) {
+// int rank;
+// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+ const size_t sz[3] = {81,81,1};
+ Box<3, double> box({0, 0,-5}, {0.5, 0.5,5});
+ size_t bc[3] = {NON_PERIODIC, NON_PERIODIC,NON_PERIODIC};
+ double spacing = box.getHigh(0) / (sz[0] - 1);
+ Ghost<3, double> ghost(spacing * 3.1);
+ double rCut = 3.1 * spacing;
+ BOOST_TEST_MESSAGE("Init vector_dist...");
+
+ vector_dist<3, double, aggregate<double,double,double,double,double,double,double[3]>> domain(0, box, bc, ghost);
+
+
//Init_DCPSE(domain)
BOOST_TEST_MESSAGE("Init domain...");
@@ -832,10 +1076,12 @@ BOOST_AUTO_TEST_CASE(dcpse_surface_sphere_old) {
}
domain.ghost_get<1,2,3>();
- SurfaceDerivative_xx<6> Dxx(domain, 2, rCut,6.0,support_options::ADAPTIVE_SURFACE);
- SurfaceDerivative_yy<6> Dyy(domain, 2, rCut,6.0,support_options::ADAPTIVE_SURFACE);
- SurfaceDerivative_zz<6> Dzz(domain, 2, rCut,6.0,support_options::ADAPTIVE_SURFACE);
-
+ SurfaceDerivative_xx<6> Dxx(domain, 2, rCut,1.0,support_options::LOAD);
+ SurfaceDerivative_yy<6> Dyy(domain, 2, rCut,1.0,support_options::LOAD);
+ SurfaceDerivative_zz<6> Dzz(domain, 2, rCut,1.0,support_options::LOAD);
+ Dxx.load(domain,"Sdxx_test");
+ Dyy.load(domain,"Sdyy_test");
+ Dzz.load(domain,"Sdzz_test");
domain.ghost_get<2>();
sol=Dxx(anasol)+Dyy(anasol)+Dzz(anasol);
@@ -854,10 +1100,12 @@ BOOST_AUTO_TEST_CASE(dcpse_surface_sphere_old) {
//std::cout << "Maximum Analytic Error: " << worst1 << std::endl;
//domain.ghost_get<4>();
- domain.write("Robin_anasol");
+ //domain.write("Robin_anasol");
BOOST_REQUIRE(worst1 < 0.03);
}
+
+
BOOST_AUTO_TEST_SUITE_END()
#endif
diff --git a/src/DCPSE/Dcpse.hpp b/src/DCPSE/Dcpse.hpp
index bc30fb3a05803e8a70de415b4441eacf7bfb5b4d..c7ca10234ac2d1b1225a0d72d65bdfb600b766a9 100644
--- a/src/DCPSE/Dcpse.hpp
+++ b/src/DCPSE/Dcpse.hpp
@@ -213,7 +213,6 @@ public:
particles.ghost_get_subset(); // This communicates which ghost particles to be excluded from support
if(opt==support_options::ADAPTIVE_SURFACE) {
- supportSizeFactor=nSpacing;
if(dim==2){
nCount=3;
}
@@ -226,18 +225,22 @@ public:
auto it = particlesTo.getDomainAndGhostIterator();
while (it.isNext()) {
auto key_o = particlesTo.getOriginKey(it.get());
- Support support = supportBuilder.getSupport(it,1,opt);
+ Support support = supportBuilder.getSupport(it,nSpacing,opt);
nSpacings.add(supportBuilder.getLastAvgspacing());
++it;
}
}
- createNormalParticles<NORMAL_ID>(particles);
+ if(opt!=support_options::LOAD) {
+ createNormalParticles<NORMAL_ID>(particles);
#ifdef SE_CLASS1
- particles.write("WithNormalParticlesQC");
+ particles.write("WithNormalParticlesQC");
#endif
+ }
initializeStaticSize(particles, particles, convergenceOrder, rCut, supportSizeFactor);
- accumulateAndDeleteNormalParticles(particles);
+ if(opt!=support_options::LOAD) {
+ accumulateAndDeleteNormalParticles(particles);
+ }
}
Dcpse(vector_type &particles,
@@ -728,61 +731,6 @@ public:
}
- /**
- * Computes the value of the Surface differential operator for one particle for o1 representing a scalar
- *
- * \param key particle
- * \param o1 source property
- * \return the selected derivative
- *
- */
- template<typename op_type>
- auto computeSurfaceDifferentialOperator(const vect_dist_key_dx &key,
- op_type &o1) -> decltype(is_scalar<std::is_fundamental<decltype(o1.value(
- key))>::value>::analyze(key, o1)) {
-
- typedef decltype(is_scalar<std::is_fundamental<decltype(o1.value(key))>::value>::analyze(key, o1)) expr_type;
-
- T sign = 1.0;
- if (differentialOrder % 2 == 0) {
- sign = -1;
- }
-
- double eps = localEps.get(key.getKey());
- double epsInvPow = localEpsInvPow.get(key.getKey());
-
- auto &particles = o1.getVector();
-
-#ifdef SE_CLASS1
- if(particles.getMapCtr()!=this->getUpdateCtr())
- {
- std::cerr<<__FILE__<<":"<<__LINE__<<" Error: You forgot a DCPSE operator update after map."<<std::endl;
- }
-#endif
-
- expr_type Dfxp = 0;
- Support support = localSupports.get(key.getKey());
- size_t xpK = support.getReferencePointKey();
- //Point<dim, T> xp = particles.getPos(xpK);
- expr_type fxp = sign * o1.value(key);
- size_t kerOff = kerOffsets.get(xpK);
- auto & keys = support.getKeys();
- for (int i = 0 ; i < keys.size() ; i++)
- {
- size_t xqK = keys.get(i);
- expr_type fxq = o1.value(vect_dist_key_dx(xqK));
- Dfxp = Dfxp + (fxq + fxp) * calcKernels.get(kerOff+i);
- }
- //additional contribution of particles normal to reference Particle
- Dfxp = Dfxp + (o1.value(key)+fxp) * calcKernels.get(kerOff+keys.size());
- Dfxp = Dfxp * epsInvPow;
- //
- //T trueDfxp = particles.template getProp<2>(xpK);
- // Store Dfxp in the right position
- return Dfxp;
- }
-
-
void initializeUpdate(vector_type &particlesFrom,vector_type2 &particlesTo)
{
#ifdef SE_CLASS1
diff --git a/src/DCPSE/SupportBuilder.hpp b/src/DCPSE/SupportBuilder.hpp
index 10e75fb397c5ff5fba7dac3bb506a9b9bb8a8bd5..404ea3afa566ca13e693b73da61fdd04af4dd418 100644
--- a/src/DCPSE/SupportBuilder.hpp
+++ b/src/DCPSE/SupportBuilder.hpp
@@ -24,54 +24,56 @@ enum support_options
template<typename vector_type,typename vector_type2>
-class SupportBuilder
-{
+class SupportBuilder {
private:
vector_type &domainFrom;
vector_type2 &domainTo;
decltype(std::declval<vector_type>().getCellList(0.0)) cellList;
const Point<vector_type::dims, unsigned int> differentialSignature;
- typename vector_type::stype rCut,AvgSpacing;
+ typename vector_type::stype rCut, AvgSpacing;
bool is_interpolation;
public:
- SupportBuilder(vector_type &domainFrom,vector_type2 &domainTo, const Point<vector_type::dims, unsigned int> differentialSignature,
- typename vector_type::stype rCut,
- bool is_interpolation)
- :domainFrom(domainFrom),
- domainTo(domainTo),
- differentialSignature(differentialSignature),
- rCut(rCut),is_interpolation(is_interpolation)
- {
+ SupportBuilder(vector_type &domainFrom, vector_type2 &domainTo,
+ const Point<vector_type::dims, unsigned int> differentialSignature,
+ typename vector_type::stype rCut,
+ bool is_interpolation)
+ : domainFrom(domainFrom),
+ domainTo(domainTo),
+ differentialSignature(differentialSignature),
+ rCut(rCut), is_interpolation(is_interpolation) {
cellList = domainFrom.getCellList(rCut);
}
- SupportBuilder(vector_type &domainFrom,vector_type2 &domainTo, unsigned int differentialSignature[vector_type::dims], typename vector_type::stype rCut, bool is_interpolation)
- : SupportBuilder(domainFrom,domainTo, Point<vector_type::dims, unsigned int>(differentialSignature), rCut) {}
+ SupportBuilder(vector_type &domainFrom, vector_type2 &domainTo,
+ unsigned int differentialSignature[vector_type::dims], typename vector_type::stype rCut,
+ bool is_interpolation)
+ : SupportBuilder(domainFrom, domainTo, Point<vector_type::dims, unsigned int>(differentialSignature),
+ rCut) {}
template<typename iterator_type>
- Support getSupport(iterator_type itPoint, unsigned int requiredSize, support_options opt)
- {
+ Support getSupport(iterator_type itPoint, unsigned int requiredSize, support_options opt) {
// Get spatial position from point iterator
vect_dist_key_dx p = itPoint.get();
vect_dist_key_dx pOrig = itPoint.getOrig();
Point<vector_type::dims, typename vector_type::stype> pos = domainTo.getPos(p.getKey());
// Get cell containing current point and add it to the set of cell keys
- grid_key_dx<vector_type::dims> curCellKey = cellList.getCellGrid(pos); // Here get the key of the cell where the current point is
+ grid_key_dx<vector_type::dims> curCellKey = cellList.getCellGrid(
+ pos); // Here get the key of the cell where the current point is
std::set<grid_key_dx<vector_type::dims>> supportCells;
supportCells.insert(curCellKey);
// Make sure to consider a set of cells providing enough points for the support
- enlargeSetOfCellsUntilSize(supportCells, requiredSize + 1,opt); // NOTE: this +1 is because we then remove the point itself
+ enlargeSetOfCellsUntilSize(supportCells, requiredSize + 1,
+ opt); // NOTE: this +1 is because we then remove the point itself
// Now return all the points from the support into a vector
- std::vector<size_t> supportKeys = getPointsInSetOfCells(supportCells,p,pOrig,requiredSize,opt);
+ std::vector<size_t> supportKeys = getPointsInSetOfCells(supportCells, p, pOrig, requiredSize, opt);
- if (is_interpolation == false)
- {
+ if (is_interpolation == false) {
auto p_o = domainFrom.getOriginKey(p.getKey());
std::remove(supportKeys.begin(), supportKeys.end(), p_o.getKey());
}
@@ -80,78 +82,65 @@ public:
return Support(p_o.getKey(), openfpm::vector_std<size_t>(supportKeys.begin(), supportKeys.end()));
}
- typename vector_type::stype getLastAvgspacing()
- {
+ typename vector_type::stype getLastAvgspacing() {
return this->AvgSpacing;
}
private:
- size_t getCellLinId(const grid_key_dx<vector_type::dims> &cellKey)
- {
+ size_t getCellLinId(const grid_key_dx<vector_type::dims> &cellKey) {
mem_id id = cellList.getGrid().LinId(cellKey);
return static_cast<size_t>(id);
}
- size_t getNumElementsInCell(const grid_key_dx<vector_type::dims> &cellKey)
- {
+ size_t getNumElementsInCell(const grid_key_dx<vector_type::dims> &cellKey) {
const size_t curCellId = getCellLinId(cellKey);
size_t numElements = cellList.getNelements(curCellId);
return numElements;
}
- size_t getNumElementsInSetOfCells(const std::set<grid_key_dx<vector_type::dims>> &set)
- {
+ size_t getNumElementsInSetOfCells(const std::set<grid_key_dx<vector_type::dims>> &set) {
size_t tot = 0;
- for (const auto cell : set)
- {
+ for (const auto cell: set) {
tot += getNumElementsInCell(cell);
}
return tot;
}
void enlargeSetOfCellsUntilSize(std::set<grid_key_dx<vector_type::dims>> &set, unsigned int requiredSize,
- support_options opt)
- {
- if (opt==support_options::RADIUS){
- auto cell=*set.begin();
+ support_options opt) {
+ if (opt == support_options::RADIUS || opt == support_options::ADAPTIVE_SURFACE) {
+ auto cell = *set.begin();
grid_key_dx<vector_type::dims> middle;
- int n=std::ceil(rCut/cellList.getCellBox().getHigh(0));
+ int n = std::ceil(rCut / cellList.getCellBox().getHigh(0));
size_t sz[vector_type::dims];
- for (int i=0;i<vector_type::dims;i++)
- {
- sz[i]=2*n+1;
- middle.set_d(i,n);
+ for (int i = 0; i < vector_type::dims; i++) {
+ sz[i] = 2 * n + 1;
+ middle.set_d(i, n);
}
- grid_sm<vector_type::dims,void> g(sz);
+ grid_sm<vector_type::dims, void> g(sz);
grid_key_dx_iterator<vector_type::dims> g_k(g);
- while(g_k.isNext())
- {
- auto key=g_k.get();
- key=cell+key-middle;
- if (isCellKeyInBounds(key))
- {
+ while (g_k.isNext()) {
+ auto key = g_k.get();
+ key = cell + key - middle;
+ if (isCellKeyInBounds(key)) {
set.insert(key);
}
++g_k;
}
- }
- else{
- while (getNumElementsInSetOfCells(set) < 5.0*requiredSize) //Why 5*requiredSize? Becasue it can help with adaptive resolutions.
+ } else {
+ while (getNumElementsInSetOfCells(set) <
+ 5.0 * requiredSize) //Why 5*requiredSize? Becasue it can help with adaptive resolutions.
{
auto tmpSet = set;
- for (const auto el : tmpSet)
- {
- for (unsigned int i = 0; i < vector_type::dims; ++i)
- {
+ for (const auto el: tmpSet) {
+ for (unsigned int i = 0; i < vector_type::dims; ++i) {
const auto pOneEl = el.move(i, +1);
const auto mOneEl = el.move(i, -1);
- if (isCellKeyInBounds(pOneEl))
- {
+ if (isCellKeyInBounds(pOneEl)) {
set.insert(pOneEl);
}
- if (isCellKeyInBounds(mOneEl))
- {
+ if (isCellKeyInBounds(mOneEl)) {
set.insert(mOneEl);
}
}
@@ -162,77 +151,73 @@ private:
}
std::vector<size_t> getPointsInSetOfCells(std::set<grid_key_dx<vector_type::dims>> set,
- vect_dist_key_dx & p,
- vect_dist_key_dx & pOrig,
- size_t requiredSupportSize,
- support_options opt)
- {
- struct reord
- {
+ vect_dist_key_dx &p,
+ vect_dist_key_dx &pOrig,
+ size_t requiredSupportSize,
+ support_options opt) {
+ struct reord {
typename vector_type::stype dist;
size_t offset;
- bool operator<(const reord & p) const
- {return this->dist < p.dist;}
+ bool operator<(const reord &p) const { return this->dist < p.dist; }
};
openfpm::vector<reord> rp;
std::vector<size_t> points;
- Point<vector_type::dims,typename vector_type::stype> xp = domainTo.getPos(p);
- for (const auto cellKey : set)
- {
+ Point<vector_type::dims, typename vector_type::stype> xp = domainTo.getPos(p);
+ for (const auto cellKey: set) {
const size_t cellLinId = getCellLinId(cellKey);
const size_t elemsInCell = getNumElementsInCell(cellKey);
- for (size_t k = 0; k < elemsInCell; ++k)
- {
+ for (size_t k = 0; k < elemsInCell; ++k) {
size_t el = cellList.get(cellLinId, k);
- if (pOrig.getKey() == el && is_interpolation == false) {continue;}
+ if (pOrig.getKey() == el && is_interpolation == false) { continue; }
- Point<vector_type::dims,typename vector_type::stype> xq = domainFrom.getPosOrig(el);
+ Point<vector_type::dims, typename vector_type::stype> xq = domainFrom.getPosOrig(el);
//points.push_back(el);
reord pr;
pr.dist = xp.distance(xq);
pr.offset = el;
-
rp.add(pr);
}
}
- if (opt == support_options::RADIUS)
- {
- for (int i = 0 ; i < rp.size() ; i++)
- {
- if (rp.get(i).dist < rCut)
- {
+ if (opt == support_options::RADIUS) {
+ for (int i = 0; i < rp.size(); i++) {
+ if (rp.get(i).dist < rCut) {
points.push_back(rp.get(i).offset);
}
}
- /* #ifdef SE_CLASS1
- if (points.size()<requiredSupportSize)
- {
- std::cerr<<__FILE__<<":"<<__LINE__<<"Note that the DCPSE neighbourhood doesn't have asked no. particles (Increase the rCut or reduce the over_sampling factor)";
- std::cout<<"Particels asked (minimum*oversampling_factor): "<<requiredSupportSize<<". Particles Possible with given options:"<<points.size()<<"."<<std::endl;
- }
- #endif*/
+ /* #ifdef SE_CLASS1
+ if (points.size()<requiredSupportSize)
+ {
+ std::cerr<<__FILE__<<":"<<__LINE__<<"Note that the DCPSE neighbourhood doesn't have asked no. particles (Increase the rCut or reduce the over_sampling factor)";
+ std::cout<<"Particels asked (minimum*oversampling_factor): "<<requiredSupportSize<<". Particles Possible with given options:"<<points.size()<<"."<<std::endl;
+ }
+ #endif*/
}
- else
- { rp.sort();
- AvgSpacing=rp.get(1).dist;
- for (int i = 0 ; i < requiredSupportSize ; i++)
- {
- if(opt==support_options::ADAPTIVE_SURFACE && rp.get(i).dist > rCut)
- {}
- else{
- points.push_back(rp.get(i).offset);
+ else if(opt == support_options::ADAPTIVE_SURFACE) {
+ AvgSpacing = std::numeric_limits<double>::max();
+ for (int i = 0; i < rp.size(); i++) {
+ if (AvgSpacing > rp.get(i).dist && rp.get(i).dist != 0) {
+ AvgSpacing = rp.get(i).dist;
+ }
+ }
+ for (int i = 0; i < rp.size(); i++) {
+ if (rp.get(i).dist < 3.9 * AvgSpacing) {
+ points.push_back(rp.get(i).offset);
}
- //AvgSpacing+=rp.get(i).dist;
}
- //AvgSpacing=AvgSpacing/requiredSupportSize;
}
-
+ else {
+ rp.sort();
+ for (int i = 0; i < requiredSupportSize; i++) {
+ points.push_back(rp.get(i).offset);
+ }
+ }
+ //AvgSpacing=AvgSpacing/requiredSupportSize
return points;
}