Commit 06ddaea0 authored by Nathan Coppin's avatar Nathan Coppin
Browse files

solving problem

parent 3bc018a9
Pipeline #8250 failed with stages
in 3 minutes and 30 seconds
......@@ -333,7 +333,6 @@ class ParticleProblem :
nameb = name.encode()
fdata.append((nameb,v[:,[0]]))
fdata.append((nameb+b"_t",v[:,[1]]))
fdata.append((nameb+b"_dir_n",basis[:,0,:self._dim-1]))
if self._dim == 3:
fdata.append((nameb+b"_s",v[:,[2]]))
fdata.append((nameb+b"_idx",o))
......@@ -381,25 +380,10 @@ class ParticleProblem :
v = []
v.append(np.vstack([fdata[k] for k in ks]))
v.append(np.vstack([fdata[k+"_t"] for k in ks]))
basis = []
cheat = np.vstack([fdata[k+"_dir_n"] for k in ks])[:,:2]
if self._dim == 2:
cheat = np.append(cheat,((1 - cheat[:,0]**2)**0.5).reshape((cheat.shape[0],1)),axis=1)
cheat[np.isnan(cheat[:,1]),1] = 0
basis.append(cheat)
basis.append(np.vstack([fdata[k+"_dir_n"] for k in ks])[:,[1,0]]*np.array([-1, 1]))
else:
cheat = np.append(cheat,((1-cheat[:,0]**2-cheat[:,1]**2)**0.5).reshape((cheat.shape[0],1)),axis=1)
cheat[np.isnan(cheat[:,2]),2] = 0
basis.append(cheat)
dirt = cheat[:,[2,0,1]]*np.array([-1,1,-1])-cheat*np.einsum("ij,ij->i",cheat,cheat[:,[2,0,1]]*np.array([-1, 1, -1]))[:,np.newaxis]
nnorm = np.linalg.norm(dirt,axis=1)
nnorm[nnorm==0] = 1
basis.append(dirt/nnorm[:,np.newaxis])
basis.append(np.cross(cheat,dirt/nnorm[:,np.newaxis]))
if self._dim == 3:
v.append(np.vstack([fdata[k+"_s"] for k in ks]))
basis = np.column_stack(basis)
v = np.column_stack(v)*contact_factor
basis = np.zeros((v.shape[0],self._dim,self._dim))
o = np.vstack([fdata[k+"_idx"] for k in ks])
t = np.repeat([0,1,2,3],[fdata[k].shape[0] for k in ks])
self._lib.particleProblemSetContacts(self._p,c_int(t.shape[0]),_np2c(o,np.uint64),_np2c(v),_np2c(basis),_np2c(t,np.int32))
......
......@@ -1121,13 +1121,6 @@ int particleProblemSolve(ParticleProblem *p, double alert, double dt, double tol
int particleProblemIterate(ParticleProblem *p, double alert, double dt, double tol, int maxit,int force_motion, double *externalForces)
{
for (size_t j = 0; j<vectorSize(p->contacts); ++j){
p->contacts[j].dv *= dt;
p->contacts[j].ct *= dt;
#if DIMENSION==3
p->contacts[j].cs *= dt;
#endif
}
double *oldVelocity = malloc(sizeof(double)*vectorSize(p->velocity));
for (size_t j = 0; j < vectorSize(p->particles); ++j){
for (size_t i = 0; i < DIMENSION; ++i){
......@@ -1149,13 +1142,6 @@ int particleProblemIterate(ParticleProblem *p, double alert, double dt, double t
}
}
free(oldVelocity);
for (size_t j = 0; j<vectorSize(p->contacts); ++j){
p->contacts[j].dv /= dt;
p->contacts[j].ct /= dt;
#if DIMENSION==3
p->contacts[j].cs /= dt;
#endif
}
return 0;
}
for (size_t i = 0; i < vectorSize(p->position); ++i){
......@@ -1183,13 +1169,6 @@ int particleProblemIterate(ParticleProblem *p, double alert, double dt, double t
p->contactForces[j*DIMENSION+i] = (p->velocity[j*DIMENSION+i]-p->contactForces[j*DIMENSION+i])*p->particles[j].m/dt;
}
}
for (size_t j = 0; j<vectorSize(p->contacts); ++j){
p->contacts[j].dv /= dt;
p->contacts[j].ct /= dt;
#if DIMENSION==3
p->contacts[j].cs /= dt;
#endif
}
free(oldVelocity);
return 1;
}
......
Markdown is supported
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