@@ -34,7 +34,7 @@ AreaConservationMeshForceCompute::AreaConservationMeshForceCompute(
34
34
if (m_ignore_type)
35
35
n_types = 1 ;
36
36
37
- GPUArray<Scalar2 > params (n_types, m_exec_conf);
37
+ GPUArray<area_conservation_param_t > params (n_types, m_exec_conf);
38
38
m_params.swap (params);
39
39
40
40
GPUArray<Scalar> area (n_types, m_exec_conf);
@@ -47,32 +47,31 @@ AreaConservationMeshForceCompute::~AreaConservationMeshForceCompute()
47
47
}
48
48
49
49
/* ! \param type Type of the angle to set parameters for
50
- \param K Stiffness parameter for the force computation
51
- \param A0 desired surface area to maintain for the force computation
50
+ \param params Parameters to set
52
51
53
52
Sets parameters for the potential of a particular mesh type
54
53
*/
55
- void AreaConservationMeshForceCompute::setParams (unsigned int type, Scalar K, Scalar A0)
54
+ void AreaConservationMeshForceCompute::setParams (unsigned int type,
55
+ const area_conservation_param_t & params)
56
56
{
57
57
if (!m_ignore_type || type == 0 )
58
58
{
59
- ArrayHandle<Scalar2> h_params (m_params, access_location::host, access_mode::readwrite);
60
- // update the local copy of the memory
61
- h_params.data [type] = make_scalar2 (K, A0);
59
+ ArrayHandle<area_conservation_param_t > h_params (m_params,
60
+ access_location::host,
61
+ access_mode::readwrite);
62
+ h_params.data [type] = params;
62
63
63
- // check for some silly errors a user could make
64
- if (K <= 0 )
64
+ if (params.k <= 0 )
65
65
m_exec_conf->msg ->warning () << " area: specified K <= 0" << endl;
66
- if (A0 <= 0 )
66
+ if (params. A0 <= 0 )
67
67
m_exec_conf->msg ->warning () << " area: specified A0 <= 0" << endl;
68
68
}
69
69
}
70
70
71
71
void AreaConservationMeshForceCompute::setParamsPython (std::string type, pybind11::dict params)
72
72
{
73
73
auto typ = m_mesh_data->getMeshBondData ()->getTypeByName (type);
74
- auto _params = area_conservation_params (params);
75
- setParams (typ, _params.k , _params.A0 );
74
+ setParams (typ, area_conservation_param_t (params));
76
75
}
77
76
78
77
pybind11::dict AreaConservationMeshForceCompute::getParams (std::string type)
@@ -85,11 +84,10 @@ pybind11::dict AreaConservationMeshForceCompute::getParams(std::string type)
85
84
}
86
85
if (m_ignore_type)
87
86
typ = 0 ;
88
- ArrayHandle<Scalar2> h_params (m_params, access_location::host, access_mode::read);
89
- pybind11::dict params;
90
- params[" k" ] = h_params.data [typ].x ;
91
- params[" A0" ] = h_params.data [typ].y ;
92
- return params;
87
+ ArrayHandle<area_conservation_param_t > h_params (m_params,
88
+ access_location::host,
89
+ access_mode::read);
90
+ return h_params.data [typ].asDict ();
93
91
}
94
92
95
93
/* ! Actually perform the force computation
@@ -107,7 +105,9 @@ void AreaConservationMeshForceCompute::computeForces(uint64_t timestep)
107
105
ArrayHandle<Scalar4> h_force (m_force, access_location::host, access_mode::overwrite);
108
106
ArrayHandle<Scalar> h_virial (m_virial, access_location::host, access_mode::overwrite);
109
107
size_t virial_pitch = m_virial.getPitch ();
110
- ArrayHandle<Scalar2> h_params (m_params, access_location::host, access_mode::read);
108
+ ArrayHandle<area_conservation_param_t > h_params (m_params,
109
+ access_location::host,
110
+ access_mode::read);
111
111
ArrayHandle<Scalar> h_area (m_area, access_location::host, access_mode::read);
112
112
113
113
ArrayHandle<typename Angle::members_t > h_triangles (
@@ -205,12 +205,13 @@ void AreaConservationMeshForceCompute::computeForces(uint64_t timestep)
205
205
else
206
206
triN = h_pts.data [triangle_type];
207
207
208
- Scalar AreaDiff = h_area.data [triangle_type] - h_params.data [triangle_type].y ;
208
+ Scalar AreaDiff = h_area.data [triangle_type] - h_params.data [triangle_type].A0 ;
209
209
210
- Scalar energy = h_params.data [triangle_type].x * AreaDiff * AreaDiff
211
- / (6 * h_params.data [triangle_type].y * triN);
210
+ Scalar energy = h_params.data [triangle_type].k * AreaDiff * AreaDiff
211
+ / (6 * h_params.data [triangle_type].A0 * triN);
212
212
213
- AreaDiff = h_params.data [triangle_type].x / h_params.data [triangle_type].y * AreaDiff / 2.0 ;
213
+ AreaDiff
214
+ = h_params.data [triangle_type].k / h_params.data [triangle_type].A0 * AreaDiff / 2.0 ;
214
215
215
216
Fab = AreaDiff * (-nab * rac * s_baac + ds_drab * rab * rac);
216
217
Fac = AreaDiff * (-nac * rab * s_baac + ds_drac * rab * rac);
0 commit comments