diff --git a/src/DCPSE_op/DCPSEActiveGelTest.cpp b/src/DCPSE_op/DCPSEActiveGelTest.cpp
index 896db418040b62238fcb059401d5e4b23f09c3e0..a8a57a9a1e080a4083005671d53b4c52d833db32 100644
--- a/src/DCPSE_op/DCPSEActiveGelTest.cpp
+++ b/src/DCPSE_op/DCPSEActiveGelTest.cpp
@@ -26,15 +26,15 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
double Ly = box.getHigh(1);
size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
double spacing = box.getHigh(0) / (sz[0] - 1);
- double rCut = 3.1 * spacing;
- double rCut2 = 3.1 * spacing;
+ double rCut = 3.2 * spacing;
+ double rCut2 = 3.2 * spacing;
int ord = 2;
double sampling_factor = 1.9;
int ord2 = 2;
double sampling_factor2 = 1.9;
Ghost<2, double> ghost(rCut);
/* pol V vort Ext Press strain stress Mfield, dPol dV RHS f1 f2 f3 f4 f5 f6 H V_t div H_t delmu */
- vector_dist<2, double, aggregate<VectorS<2, double>, VectorS<2, double>, double[2][2], VectorS<2, double>, double, double[2][2], double[2][2], VectorS<2, double>, VectorS<2, double>, VectorS<2, double>, VectorS<2, double>, double, double, double, double, double, double, double, VectorS<2, double>, double, double, double[2], double[2], double[2], double[2], double[2], double[2], double,VectorS<2, double>,VectorS<2, double>,VectorS<2, double>,VectorS<2, double>>> Particles(
+ vector_dist<2, double, aggregate<VectorS<2, double>, VectorS<2, double>, double[2][2], VectorS<2, double>, double, double[2][2], double[2][2], VectorS<2, double>, VectorS<2, double>, VectorS<2, double>, VectorS<2, double>, double, double, double, double, double, double, double, VectorS<2, double>, double, double, double[2], double[2], double[2], double[2], double[2], double[2], double, VectorS<2, double>, VectorS<2, double>, VectorS<2, double>, VectorS<2, double>>> Particles(
0, box, bc, ghost);
auto it = Particles.getGridIterator(sz);
@@ -70,7 +70,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
constexpr int Strain_rate = 5;
constexpr int Stress = 6;
constexpr int MolField = 7;
- auto Pos=getV<PROP_POS>(Particles);
+ auto Pos = getV<PROP_POS>(Particles);
auto Pol = getV<Polarization>(Particles);
auto V = getV<Velocity>(Particles);
@@ -192,19 +192,31 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
Laplacian Lap(Particles, ord, rCut, sampling_factor, support_options::RADIUS);
Divergence Div(Particles, ord2, rCut2, sampling_factor2, support_options::RADIUS);
-
+ petsc_solver<double> solverPetsc;
+ //solverPetsc.setSolver(KSPGMRES);
+ //solverPetsc.setRestart(250);
+ //solverPetsc.setPreconditioner(PCJACOBI);
+ petsc_solver<double> solverPetsc2;
+ //solverPetsc2.setSolver(KSPGMRES);
+ //solverPetsc2.setRestart(250);
+ /*solverPetsc2.setPreconditioner(PCJACOBI);
+ solverPetsc2.setRelTol(1e-6);
+ solverPetsc2.setAbsTol(1e-5);
+ solverPetsc2.setDivTol(1e4);*/
eq_id vx, vy;
timer tt;
+ timer tt3;
vx.setId(0);
vy.setId(1);
- int n = 50;
- double dt = 1e-3;
- double tim=0;
- double tf=1e-2;
- while (tim <= tf)
- {
- Particles.ghost_get<Polarization>();
+ int n = 75;
+ int ctr = 0;
+ double dt = 3e-3;
+ double tim = 0;
+ double tf = 7.5;
+ while (tim <= tf) {
+ tt.start();
+ //Particles.ghost_get<Polarization>();
sigma[x][x] =
-Ks * Dx(Pol[x]) * Dx(Pol[x]) - Kb * Dx(Pol[y]) * Dx(Pol[y]) + (Kb - Ks) * Dy(Pol[x]) * Dx(Pol[y]);
sigma[x][y] =
@@ -218,7 +230,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
h[y] = (Pol[x] * (Ks * Dyy(Pol[y]) + Kb * Dxx(Pol[y]) + (Ks - Kb) * Dxy(Pol[x])) -
Pol[y] * (Ks * Dxx(Pol[x]) + Kb * Dyy(Pol[x]) + (Ks - Kb) * Dxy(Pol[y])));
- Particles.ghost_get<MolField>();
+ //Particles.ghost_get<MolField>();
f1 = gama * nu * Pol[x] * Pol[x] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) /
(Pol[x] * Pol[x] + Pol[y] * Pol[y]);
@@ -229,7 +241,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
f4 = 2.0 * gama * nu * Pol[x] * Pol[x] * Pol[x] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y]);
f5 = 4.0 * gama * nu * Pol[x] * Pol[x] * Pol[y] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y]);
f6 = 2.0 * gama * nu * Pol[x] * Pol[y] * Pol[y] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y]);
- Particles.ghost_get<11, 12, 13, 14, 15, 16>();
+ //Particles.ghost_get<11, 12, 13, 14, 15, 16>();
Df1[x] = Dx(f1);
Df2[x] = Dx(f2);
Df3[x] = Dx(f3);
@@ -243,7 +255,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
Df4[y] = Dy(f4);
Df5[y] = Dy(f5);
Df6[y] = Dy(f6);
- Particles.ghost_get<21, 22, 23, 24, 25, 26>();
+ //Particles.ghost_get<21, 22, 23, 24, 25, 26>();
dV[x] = -0.5 * Dy(h[y]) + zeta * Dx(delmu * Pol[x] * Pol[x]) + zeta * Dy(delmu * Pol[x] * Pol[y]) -
@@ -262,12 +274,12 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
g[y]
- 0.5 * nu * Dy(gama * lambda * delmu * (Pol[x] * Pol[x] - Pol[y] * Pol[y]))
- 0.5 * Dx(-2.0 * gama * lambda * delmu * (Pol[x] * Pol[y]));
- Particles.ghost_get<9>();
+ //Particles.ghost_get<9>();
- Particles.write_frame("Polar10", 0);
+ Particles.write("PolarI");
//Velocity Solution n iterations
- double sum = 0, sum1 = 0;
+ double sum, sum1;
auto Stokes1 = eta * Lap(V[x])
+ 0.5 * nu * (Df1[x] * Dx(V[x]) + f1 * Dxx(V[x]))
@@ -283,20 +295,9 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
+ 0.5 * nu * (Df4[x] * Dx(V[x]) + f4 * Dxx(V[x]))
+ 0.5 * nu * (Df5[x] * 0.5 * (Dx(V[y]) + Dy(V[x])) + f5 * 0.5 * (Dxx(V[y]) + Dyx(V[x])))
+ 0.5 * nu * (Df6[x] * Dy(V[y]) + f6 * Dyx(V[y]));
-
- petsc_solver<double> solverPetsc;
- //solverPetsc.setSolver(KSPGMRES);
- //solverPetsc.setRestart(250);
- //solverPetsc.setPreconditioner(PCJACOBI);
- petsc_solver<double> solverPetsc2;
- //solverPetsc2.setSolver(KSPGMRES);
- //solverPetsc2.setRestart(250);
- /*solverPetsc2.setPreconditioner(PCJACOBI);
- solverPetsc2.setRelTol(1e-6);
- solverPetsc2.setAbsTol(1e-5);
- solverPetsc2.setDivTol(1e4);*/
+ tt.stop();
+ std::cout << "Init of Velocity took " << tt.getwct() << " seconds." << std::endl;
tt.start();
-
for (int i = 1; i <= n; i++) {
RHS[x] = Dx(P) + dV[x];
RHS[y] = Dy(P) + dV[y];
@@ -326,9 +327,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
SolverH.impose(H, dw_p, 0);
SolverH.impose(H, l_p, 0);
SolverH.impose(H, r_p, 0);
- tt.start();
SolverH.solve_with_solver(solverPetsc2, H);
- tt.stop();
//std::cout << "Helmholtz Solved in " << tt.getwct() << " seconds." << std::endl;
Particles.ghost_get<17>();
Particles.ghost_get<Velocity>();
@@ -390,50 +389,89 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
W[y][y] = 0;
- Particles.write_frame("Polar_Petsc", int(tim / dt));
+ Particles.write_frame("Polar_Petsc", ctr);
+ ctr++;
+
h[x] = -gama * lambda * delmu + u[x][x] * gama * nu * Pol[x] * Pol[x] / (Pol[x] * Pol[x] + Pol[y] * Pol[y])
+ u[y][y] * gama * nu * Pol[y] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y])
+ 2 * u[x][y] * gama * nu * Pol[x] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y]);
- Pos=Pos+dt*V;
+ Pos = Pos + dt * V;
tt.start();
- Derivative_x Dx(Particles, ord2, rCut2, sampling_factor2, support_options::RADIUS);
- Derivative_y Dy(Particles, ord2, rCut2, sampling_factor2, support_options::RADIUS);
- Derivative_xy Dxy(Particles, ord, rCut, sampling_factor, support_options::RADIUS);
+ Dx.update(Particles);
+ Dy.update(Particles);
+ Dxy.update(Particles);
auto Dyx = Dxy;
- Derivative_xx Dxx(Particles, ord, rCut, sampling_factor, support_options::RADIUS);
- Derivative_yy Dyy(Particles, ord, rCut, sampling_factor, support_options::RADIUS);
- Laplacian Lap(Particles, ord, rCut, sampling_factor, support_options::RADIUS);
- Divergence Div(Particles, ord2, rCut2, sampling_factor2, support_options::RADIUS);
+ Dxx.update(Particles);
+ Dyy.update(Particles);
+ Lap.update(Particles);
+ Div.update(Particles);
tt.stop();
- std::cout << "Updattion of operators took " << tt.getwct() << " seconds."<< std::endl;
-
- k1[x]=((h[x] * Pol[x] - h[y] * Pol[y]) / gama + lambda * delmu * Pol[x] -
- nu * (u[x][x] * Pol[x] + u[x][y] * Pol[y]) + W[x][x] * Pol[x] + W[x][y] * Pol[y]);
- k1[y]=((h[x] * Pol[y] - h[y] * Pol[x]) / gama + lambda * delmu * Pol[y] -
- nu * (u[y][x] * Pol[x] + u[y][y] * Pol[y]) + W[y][x] * Pol[x] + W[y][y] * Pol[y]);
-
- k2[x]=((h[x] * 0.5*dt*k1[x]*Pol[x] - h[y] * 0.5*dt*k1[y]*Pol[y]) / gama + lambda * delmu * 0.5*dt*k1[x]*Pol[x] -
- nu * (u[x][x] * 0.5*dt*k1[x]*Pol[x] + u[x][y] * 0.5*dt*k1[y]*Pol[y]) + W[x][x] * 0.5*dt*k1[x]*Pol[x] + W[x][y] * 0.5*dt*k1[y]*Pol[y]);
- k2[y]=((h[x] * 0.5*dt*k1[y]*Pol[y] - h[y] * 0.5*dt*k1[x]*Pol[x]) / gama + lambda * delmu *0.5*dt*k1[y]* Pol[y] -
- nu * (u[y][x] *0.5*dt*k1[x]* Pol[x] + u[y][y] *0.5*dt*k1[y]* Pol[y]) + W[y][x] *0.5*dt*k1[x]* Pol[x] + W[y][y] *0.5*dt*k1[y]* Pol[y]);
-
- k3[x]=((h[x] * 0.5*dt*k2[x]*Pol[x] - h[y] * 0.5*dt*k2[y]*Pol[y]) / gama + lambda * delmu * 0.5*dt*k2[x]*Pol[x] -
- nu * (u[x][x] * 0.5*dt*k2[x]*Pol[x] + u[x][y] * 0.5*dt*k2[y]*Pol[y]) + W[x][x] * 0.5*dt*k2[x]*Pol[x] + W[x][y] * 0.5*dt*k2[y]*Pol[y]);
- k3[y]=((h[x] * 0.5*dt*k2[y]*Pol[y] - h[y] *0.5*dt*k2[x]* Pol[x]) / gama + lambda * delmu * 0.5*dt*k2[y]*Pol[y] -
- nu * (u[y][x] * 0.5*dt*k2[x]*Pol[x] + u[y][y] * 0.5*dt*k2[y]*Pol[y]) + W[y][x] *0.5*dt*k2[x]* Pol[x] + W[y][y] *0.5*dt*k2[y]* Pol[y]);
-
- k4[x]=((h[x] * dt*k3[x]*Pol[x] - h[y] * dt*k3[y]*Pol[y]) / gama + lambda * delmu * dt*k3[x]*Pol[x] -
- nu * (u[x][x] * dt*k3[x]*Pol[x] + u[x][y] * dt*k3[y]*Pol[y]) + W[x][x] * dt*k3[x]*Pol[x] + W[x][y] * dt*k3[y]*Pol[y]);
- k4[y]=((h[x] * dt*k3[y]*Pol[y] - h[y] * dt*k3[x]*Pol[x]) / gama + lambda * delmu * dt*k3[y]*Pol[y] -
- nu * (u[y][x] * dt*k3[x]*Pol[x] + u[y][y] * dt*k3[y]*Pol[y]) + W[y][x] * dt*k3[x]*Pol[x] + W[y][y] *dt*k3[y]* Pol[y]);
-
- Pol= Pol + (dt/6.0)*(k1+(2.0*k2)+(2.0*k3)+k4);
-
- /*Pol[x] = Pol[x] + dt * ((h[x] * Pol[x] - h[y] * Pol[y]) / gama + lambda * delmu * Pol[x] -
- nu * (u[x][x] * Pol[x] + u[x][y] * Pol[y]) + W[x][x] * Pol[x] + W[x][y] * Pol[y]);
- Pol[y] = Pol[y] + dt * ((h[x] * Pol[y] - h[y] * Pol[x]) / gama + lambda * delmu * Pol[y] -
- nu * (u[y][x] * Pol[x] + u[y][y] * Pol[y]) + W[y][x] * Pol[x] + W[y][y] * Pol[y]);*/
+ std::cout << "Updattion of operators took " << tt.getwct() << " seconds." << std::endl;
+
+ k1[x] = ((h[x] * Pol[x] - h[y] * Pol[y]) / gama + lambda * delmu * Pol[x] -
+ nu * (u[x][x] * Pol[x] + u[x][y] * Pol[y]) + W[x][x] * Pol[x] + W[x][y] * Pol[y]);
+ k1[y] = ((h[x] * Pol[y] - h[y] * Pol[x]) / gama + lambda * delmu * Pol[y] -
+ nu * (u[y][x] * Pol[x] + u[y][y] * Pol[y]) + W[y][x] * Pol[x] + W[y][y] * Pol[y]);
+
+ k2[x] = ((h[x] * 0.5 * dt * k1[x] * Pol[x] - h[y] * 0.5 * dt * k1[y] * Pol[y]) / gama +
+ lambda * delmu * 0.5 * dt * k1[x] * Pol[x] -
+ nu * (u[x][x] * 0.5 * dt * k1[x] * Pol[x] + u[x][y] * 0.5 * dt * k1[y] * Pol[y]) +
+ W[x][x] * 0.5 * dt * k1[x] * Pol[x] + W[x][y] * 0.5 * dt * k1[y] * Pol[y]);
+ k2[y] = ((h[x] * 0.5 * dt * k1[y] * Pol[y] - h[y] * 0.5 * dt * k1[x] * Pol[x]) / gama +
+ lambda * delmu * 0.5 * dt * k1[y] * Pol[y] -
+ nu * (u[y][x] * 0.5 * dt * k1[x] * Pol[x] + u[y][y] * 0.5 * dt * k1[y] * Pol[y]) +
+ W[y][x] * 0.5 * dt * k1[x] * Pol[x] + W[y][y] * 0.5 * dt * k1[y] * Pol[y]);
+
+ k3[x] = ((h[x] * 0.5 * dt * k2[x] * Pol[x] - h[y] * 0.5 * dt * k2[y] * Pol[y]) / gama +
+ lambda * delmu * 0.5 * dt * k2[x] * Pol[x] -
+ nu * (u[x][x] * 0.5 * dt * k2[x] * Pol[x] + u[x][y] * 0.5 * dt * k2[y] * Pol[y]) +
+ W[x][x] * 0.5 * dt * k2[x] * Pol[x] + W[x][y] * 0.5 * dt * k2[y] * Pol[y]);
+ k3[y] = ((h[x] * 0.5 * dt * k2[y] * Pol[y] - h[y] * 0.5 * dt * k2[x] * Pol[x]) / gama +
+ lambda * delmu * 0.5 * dt * k2[y] * Pol[y] -
+ nu * (u[y][x] * 0.5 * dt * k2[x] * Pol[x] + u[y][y] * 0.5 * dt * k2[y] * Pol[y]) +
+ W[y][x] * 0.5 * dt * k2[x] * Pol[x] + W[y][y] * 0.5 * dt * k2[y] * Pol[y]);
+
+ k4[x] = ((h[x] * dt * k3[x] * Pol[x] - h[y] * dt * k3[y] * Pol[y]) / gama +
+ lambda * delmu * dt * k3[x] * Pol[x] -
+ nu * (u[x][x] * dt * k3[x] * Pol[x] + u[x][y] * dt * k3[y] * Pol[y]) +
+ W[x][x] * dt * k3[x] * Pol[x] + W[x][y] * dt * k3[y] * Pol[y]);
+ k4[y] = ((h[x] * dt * k3[y] * Pol[y] - h[y] * dt * k3[x] * Pol[x]) / gama +
+ lambda * delmu * dt * k3[y] * Pol[y] -
+ nu * (u[y][x] * dt * k3[x] * Pol[x] + u[y][y] * dt * k3[y] * Pol[y]) +
+ W[y][x] * dt * k3[x] * Pol[x] + W[y][y] * dt * k3[y] * Pol[y]);
+
+ Pol = Pol + (dt / 6.0) * (k1 + (2.0 * k2) + (2.0 * k3) + k4);
+
+ for (int j = 0; j < up_p.size(); j++) {
+ auto p = up_p.get<0>(j);
+ Particles.getProp<0>(p)[x] = sin(2 * M_PI * (cos((2 * Particles.getPos(p)[x] - Lx) / Lx) -
+ sin((2 * Particles.getPos(p)[y] - Ly) / Ly)));
+ Particles.getProp<0>(p)[y] = cos(2 * M_PI * (cos((2 * Particles.getPos(p)[x] - Lx) / Lx) -
+ sin((2 * Particles.getPos(p)[y] - Ly) / Ly)));
+ }
+ for (int j = 0; j < dw_p.size(); j++) {
+ auto p = dw_p.get<0>(j);
+ Particles.getProp<0>(p)[x] = sin(2 * M_PI * (cos((2 * Particles.getPos(p)[x] - Lx) / Lx) -
+ sin((2 * Particles.getPos(p)[y] - Ly) / Ly)));
+ Particles.getProp<0>(p)[y] = cos(2 * M_PI * (cos((2 * Particles.getPos(p)[x] - Lx) / Lx) -
+ sin((2 * Particles.getPos(p)[y] - Ly) / Ly)));
+ }
+ for (int j = 0; j < l_p.size(); j++) {
+ auto p = l_p.get<0>(j);
+ Particles.getProp<0>(p)[x] = sin(2 * M_PI * (cos((2 * Particles.getPos(p)[x] - Lx) / Lx) -
+ sin((2 * Particles.getPos(p)[y] - Ly) / Ly)));
+ Particles.getProp<0>(p)[y] = cos(2 * M_PI * (cos((2 * Particles.getPos(p)[x] - Lx) / Lx) -
+ sin((2 * Particles.getPos(p)[y] - Ly) / Ly)));
+ }
+ for (int j = 0; j < r_p.size(); j++) {
+ auto p = r_p.get<0>(j);
+ Particles.getProp<0>(p)[x] = sin(2 * M_PI * (cos((2 * Particles.getPos(p)[x] - Lx) / Lx) -
+ sin((2 * Particles.getPos(p)[y] - Ly) / Ly)));
+ Particles.getProp<0>(p)[y] = cos(2 * M_PI * (cos((2 * Particles.getPos(p)[x] - Lx) / Lx) -
+ sin((2 * Particles.getPos(p)[y] - Ly) / Ly)));
+ }
+
std::cout << "Time step " << tim << " over." << std::endl;
tim += dt;
std::cout << "----------------------------------------------------------" << std::endl;
@@ -441,7 +479,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
}
Particles.deleteGhost();
tt2.stop();
- std::cout << "The simulation took" << tt2.getwct() << "Seconds.";
+ std::cout << "The simulation took " << tt2.getwct() << "Seconds.";
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
diff --git a/src/DCPSE_op/DCPSE_copy/Dcpse.hpp b/src/DCPSE_op/DCPSE_copy/Dcpse.hpp
index 49c20e2bef3d3ac7ccdced4c3cb8ac1c71fe6912..657b1c326cab01144874db5a2cd52b2511c4a384 100644
--- a/src/DCPSE_op/DCPSE_copy/Dcpse.hpp
+++ b/src/DCPSE_op/DCPSE_copy/Dcpse.hpp
@@ -1,6 +1,6 @@
//
// Created by tommaso on 29/03/19.
-//
+// Modified by Abhinav and Pietro
#ifndef OPENFPM_PDATA_DCPSE_HPP
#define OPENFPM_PDATA_DCPSE_HPP
@@ -397,6 +397,47 @@ public:
return Dfxp;
}
+ void initializeUpdate(vector_type &particles)
+ {
+ auto it = particles.getDomainIterator();
+ while (it.isNext()) {
+ // Get the points in the support of the DCPSE kernel and store the support for reuse
+ //Support<dim, T, part_type> support = supportBuilder.getSupport(it, requiredSupportSize,opt);
+ auto p = it.get();
+ Support<dim, T, part_type> support = localSupports[p.getKey()];
+ EMatrix<T, Eigen::Dynamic, Eigen::Dynamic> V(support.size(), monomialBasis.size());
+ // Vandermonde matrix computation
+ Vandermonde<dim, T, EMatrix<T, Eigen::Dynamic, Eigen::Dynamic>>
+ vandermonde(support, monomialBasis);
+ vandermonde.getMatrix(V);
+
+ T eps = vandermonde.getEps();
+
+ localSupports.push_back(support);
+ localEps.push_back(eps);
+ // Compute the diagonal matrix E
+ DcpseDiagonalScalingMatrix<dim> diagonalScalingMatrix(monomialBasis);
+ EMatrix<T, Eigen::Dynamic, Eigen::Dynamic> E(support.size(), support.size());
+ diagonalScalingMatrix.buildMatrix(E, support, eps);
+ // Compute intermediate matrix B
+ EMatrix<T, Eigen::Dynamic, Eigen::Dynamic> B = E * V;
+ // Compute matrix A
+ EMatrix<T, Eigen::Dynamic, Eigen::Dynamic> A = B.transpose() * B;
+ // Compute RHS vector b
+ DcpseRhs<dim> rhs(monomialBasis, differentialSignature);
+ EMatrix<T, Eigen::Dynamic, 1> b(monomialBasis.size(), 1);
+ rhs.template getVector<T>(b);
+ // Get the vector where to store the coefficients...
+ EMatrix<T, Eigen::Dynamic, 1> a(monomialBasis.size(), 1);
+ // ...solve the linear system...
+ a = A.colPivHouseholderQr().solve(b);
+ // ...and store the solution for later reuse
+ localCoefficients.push_back(a);
+ //
+ ++it;
+ }
+ }
+
private:
void initializeAdaptive(vector_type &particles,
@@ -511,6 +552,8 @@ private:
}
+
+
T computeKernel(Point<dim, T> x, EMatrix<T, Eigen::Dynamic, 1> a) const {
T res = 0;
unsigned int counter = 0;
diff --git a/src/DCPSE_op/DCPSE_op.hpp b/src/DCPSE_op/DCPSE_op.hpp
index c443dd0a76b12820877a6902c88c9a1e40e530d2..c72cdae8e68e0544961c62a158514c4e37b0924e 100644
--- a/src/DCPSE_op/DCPSE_op.hpp
+++ b/src/DCPSE_op/DCPSE_op.hpp
@@ -677,7 +677,6 @@ public:
}
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;
@@ -688,10 +687,19 @@ public:
template<unsigned int prp, typename particles_type>
void DrawKernel(particles_type &particles,int k)
{
- auto dcpse2 = (Dcpse<particles_type::dims,particles_type> *)dcpse;
- dcpse2->template DrawKernel<prp>(particles,k);
+ auto dcpse_temp = (Dcpse<particles_type::dims,particles_type> *)dcpse;
+ dcpse_temp->template DrawKernel<prp>(particles,k);
}
+
+ template<typename particles_type>
+ void update(particles_type &particles)
+ {
+ auto dcpse_temp = (Dcpse<particles_type::dims,particles_type> *)dcpse;
+ dcpse_temp-> initializeUpdate(particles);
+
+ }
+
};
@@ -735,6 +743,14 @@ public:
}
+ template<typename particles_type>
+ void update(particles_type &particles)
+ {
+ auto dcpse_temp = (Dcpse<particles_type::dims,particles_type> *)dcpse;
+ dcpse_temp-> initializeUpdate(particles);
+
+ }
+
};
@@ -763,6 +779,15 @@ public:
return vector_dist_expression_op<operand_type,dcpse_type,VECT_DCPSE>(arg,*(dcpse_type *)dcpse);
}
+
+ template<typename particles_type>
+ void update(particles_type &particles)
+ {
+ auto dcpse_temp = (Dcpse<particles_type::dims,particles_type> *)dcpse;
+ dcpse_temp-> initializeUpdate(particles);
+
+ }
+
};
@@ -814,6 +839,18 @@ public:
}
}
+
+ template<typename particles_type>
+ void update(particles_type &particles)
+ {
+ Dcpse<particles_type::dims, particles_type> *dcpse_ptr = (Dcpse<particles_type::dims, particles_type> *) dcpse;
+ for (int i = 0 ; i < particles_type::dims ; i++) {
+ dcpse_ptr[i]. initializeUpdate(particles);
+ }
+
+ }
+
+
};
class Curl2D
@@ -909,6 +946,17 @@ public:
}
+ template<typename particles_type>
+ void update(particles_type &particles)
+ {
+ Dcpse<particles_type::dims, particles_type> *dcpse_ptr = (Dcpse<particles_type::dims, particles_type> *) dcpse;
+ for (int i = 0 ; i < particles_type::dims ; i++) {
+ dcpse_ptr[i]. initializeUpdate(particles);
+ }
+
+ }
+
+
};
@@ -944,6 +992,17 @@ public:
return vector_dist_expression_op<operand_type,dcpse_type,VECT_DCPSE_V_DIV>(arg,*(dcpse_type(*)[operand_type::vtype::dims])dcpse);
}
+
+ template<typename particles_type>
+ void update(particles_type &particles)
+ {
+ Dcpse<particles_type::dims, particles_type> *dcpse_ptr = (Dcpse<particles_type::dims, particles_type> *) dcpse;
+ for (int i = 0 ; i < particles_type::dims ; i++) {
+ dcpse_ptr[i]. initializeUpdate(particles);
+ }
+
+ }
+
};
@@ -1005,6 +1064,18 @@ public:
}
}
+
+ template<typename particles_type>
+ void update(particles_type &particles)
+ {
+ Dcpse<particles_type::dims, particles_type> *dcpse_ptr = (Dcpse<particles_type::dims, particles_type> *) dcpse;
+ for (int i = 0 ; i < particles_type::dims ; i++) {
+ dcpse_ptr[i]. initializeUpdate(particles);
+ }
+
+ }
+
+
};
@@ -1051,6 +1122,14 @@ public:
dcpse_ptr[0].template DrawKernel<prp>(particles,k);
}
+
+ template<typename particles_type>
+ void update(particles_type &particles)
+ {
+ auto dcpse_temp = (Dcpse<particles_type::dims,particles_type> *)dcpse;
+ dcpse_temp-> initializeUpdate(particles);
+
+ }
};
@@ -1080,6 +1159,14 @@ public:
return vector_dist_expression_op<operand_type,dcpse_type,VECT_DCPSE>(arg,*(dcpse_type *)dcpse);
}
+
+ template<typename particles_type>
+ void update(particles_type &particles)
+ {
+ auto dcpse_temp = (Dcpse<particles_type::dims,particles_type> *)dcpse;
+ dcpse_temp-> initializeUpdate(particles);
+
+ }
};
@@ -1109,6 +1196,14 @@ public:
return vector_dist_expression_op<operand_type,dcpse_type,VECT_DCPSE>(arg,*(dcpse_type *)dcpse);
}
+
+ template<typename particles_type>
+ void update(particles_type &particles)
+ {
+ auto dcpse_temp = (Dcpse<particles_type::dims,particles_type> *)dcpse;
+ dcpse_temp-> initializeUpdate(particles);
+
+ }
};