Skip to content
Snippets Groups Projects
Commit b32c4faf authored by Abhinav Singh's avatar Abhinav Singh
Browse files

Merge remote-tracking branch 'origin/FD_solver' into develop 

# Conflicts:
#	src/DCPSE/SupportBuilder.hpp
parents d476fc1d 400afe6b
No related branches found
No related tags found
1 merge request!16Fd solver
Hide whitespace changes
Inline Side-by-side
src/DCPSE/DCPSE_op/DCPSE_surface_op.hpp
+ 289
− 0
View file @ b32c4faf
......@@ -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.
*
*
......
src/DCPSE/DCPSE_op/tests/DCPSE_op_Surface_tests.cpp
+ 253
− 5
View file @ b32c4faf
......@@ -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
......
src/DCPSE/Dcpse.hpp
+ 8
− 60
View file @ b32c4faf
......@@ -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
......
src/DCPSE/SupportBuilder.hpp
+ 84
− 99
View file @ b32c4faf
......@@ -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;
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment