1 /* ----------------------------------------------------------------------
2 This is the
3
4 ██╗ ██╗ ██████╗ ██████╗ ██████╗ ██╗ ██╗████████╗███████╗
5 ██║ ██║██╔════╝ ██╔════╝ ██╔════╝ ██║ ██║╚══██╔══╝██╔════╝
6 ██║ ██║██║ ███╗██║ ███╗██║ ███╗███████║ ██║ ███████╗
7 ██║ ██║██║ ██║██║ ██║██║ ██║██╔══██║ ██║ ╚════██║
8 ███████╗██║╚██████╔╝╚██████╔╝╚██████╔╝██║ ██║ ██║ ███████║
9 ╚══════╝╚═╝ ╚═════╝ ╚═════╝ ╚═════╝ ╚═╝ ╚═╝ ╚═╝ ╚══════╝®
10
11 DEM simulation engine, released by
12 DCS Computing Gmbh, Linz, Austria
13 http://www.dcs-computing.com, office@dcs-computing.com
14
15 LIGGGHTS® is part of CFDEM®project:
16 http://www.liggghts.com | http://www.cfdem.com
17
18 Core developer and main author:
19 Christoph Kloss, christoph.kloss@dcs-computing.com
20
21 LIGGGHTS® is open-source, distributed under the terms of the GNU Public
22 License, version 2 or later. It is distributed in the hope that it will
23 be useful, but WITHOUT ANY WARRANTY; without even the implied warranty
24 of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. You should have
25 received a copy of the GNU General Public License along with LIGGGHTS®.
26 If not, see http://www.gnu.org/licenses . See also top-level README
27 and LICENSE files.
28
29 LIGGGHTS® and CFDEM® are registered trade marks of DCS Computing GmbH,
30 the producer of the LIGGGHTS® software and the CFDEM®coupling software
31 See http://www.cfdem.com/terms-trademark-policy for details.
32
33 -------------------------------------------------------------------------
34 Contributing author and copyright for this file:
35 (if not contributing author is listed, this file has been contributed
36 by the core developer)
37
38 Copyright 2012- DCS Computing GmbH, Linz
39 Copyright 2009-2012 JKU Linz
40 ------------------------------------------------------------------------- */
41
42 #include "fix_heat_gran_conduction.h"
43
44 #include "atom.h"
45 #include "compute_pair_gran_local.h"
46 #include "fix_property_atom.h"
47 #include "fix_property_global.h"
48 #include "force.h"
49 #include "math_extra.h"
50 #include "properties.h"
51 #include "modify.h"
52 #include "neigh_list.h"
53 #include "pair_gran.h"
54 #include <cmath>
55 #include <algorithm>
56
57 using namespace LAMMPS_NS;
58 using namespace FixConst;
59
60 // modes for conduction contact area calaculation
61 // same as in fix_wall_gran.cpp
62
63 enum{ CONDUCTION_CONTACT_AREA_OVERLAP,
64 CONDUCTION_CONTACT_AREA_CONSTANT,
65 CONDUCTION_CONTACT_AREA_PROJECTION};
66
67 /* ---------------------------------------------------------------------- */
68
FixHeatGranCond(class LAMMPS * lmp,int narg,char ** arg)69 FixHeatGranCond::FixHeatGranCond(class LAMMPS *lmp, int narg, char **arg) :
70 FixHeatGran(lmp, narg, arg),
71 fix_conductivity_(0),
72 conductivity_(0),
73 store_contact_data_(false),
74 fix_conduction_contact_area_(0),
75 fix_n_conduction_contacts_(0),
76 fix_wall_heattransfer_coeff_(0),
77 fix_wall_temperature_(0),
78 conduction_contact_area_(0),
79 n_conduction_contacts_(0),
80 wall_heattransfer_coeff_(0),
81 wall_temp_(0),
82 area_calculation_mode_(CONDUCTION_CONTACT_AREA_OVERLAP),
83 fixed_contact_area_(0.),
84 area_correction_flag_(0),
85 deltan_ratio_(0)
86 {
87 iarg_ = 5;
88
89 bool hasargs = true;
90 while(iarg_ < narg && hasargs)
91 {
92 hasargs = false;
93 if(strcmp(arg[iarg_],"contact_area") == 0) {
94
95 if(strcmp(arg[iarg_+1],"overlap") == 0)
96 area_calculation_mode_ = CONDUCTION_CONTACT_AREA_OVERLAP;
97 else if(strcmp(arg[iarg_+1],"projection") == 0)
98 area_calculation_mode_ = CONDUCTION_CONTACT_AREA_PROJECTION;
99 else if(strcmp(arg[iarg_+1],"constant") == 0)
100 {
101 if (iarg_+3 > narg)
102 error->fix_error(FLERR,this,"not enough arguments for keyword 'contact_area constant'");
103 area_calculation_mode_ = CONDUCTION_CONTACT_AREA_CONSTANT;
104 fixed_contact_area_ = force->numeric(FLERR,arg[iarg_+2]);
105 if (fixed_contact_area_ <= 0.)
106 error->fix_error(FLERR,this,"'contact_area constant' value must be > 0");
107 iarg_++;
108 }
109 else error->fix_error(FLERR,this,"expecting 'overlap', 'projection' or 'constant' after 'contact_area'");
110 iarg_ += 2;
111 hasargs = true;
112 } else if(strcmp(arg[iarg_],"area_correction") == 0) {
113 if (iarg_+2 > narg) error->fix_error(FLERR,this,"not enough arguments for keyword 'area_correction'");
114 if(strcmp(arg[iarg_+1],"yes") == 0)
115 area_correction_flag_ = 1;
116 else if(strcmp(arg[iarg_+1],"no") == 0)
117 area_correction_flag_ = 0;
118 else error->fix_error(FLERR,this,"expecting 'yes' or 'no' after 'area_correction'");
119 iarg_ += 2;
120 hasargs = true;
121 } else if(strcmp(arg[iarg_],"store_contact_data") == 0) {
122 if (iarg_+2 > narg) error->fix_error(FLERR,this,"not enough arguments for keyword 'store_contact_data'");
123 if(strcmp(arg[iarg_+1],"yes") == 0)
124 store_contact_data_ = true;
125 else if(strcmp(arg[iarg_+1],"no") == 0)
126 store_contact_data_ = false;
127 else error->fix_error(FLERR,this,"expecting 'yes' or 'no' after 'store_contact_data'");
128 iarg_ += 2;
129 hasargs = true;
130 } else if(strcmp(style,"heat/gran/conduction") == 0)
131 error->fix_error(FLERR,this,"unknown keyword");
132 }
133
134 if(CONDUCTION_CONTACT_AREA_OVERLAP != area_calculation_mode_ && 1 == area_correction_flag_)
135 error->fix_error(FLERR,this,"can use 'area_correction' only for 'contact_area = overlap'");
136 }
137
138 /* ---------------------------------------------------------------------- */
139
~FixHeatGranCond()140 FixHeatGranCond::~FixHeatGranCond()
141 {
142
143 if (conductivity_)
144 delete []conductivity_;
145 }
146
147 /* ---------------------------------------------------------------------- */
148
post_create()149 void FixHeatGranCond::post_create()
150 {
151 FixHeatGran::post_create();
152
153 // register contact storage
154 fix_conduction_contact_area_ = static_cast<FixPropertyAtom*>(modify->find_fix_property("contactAreaConduction","property/atom","scalar",0,0,this->style,false));
155 if(!fix_conduction_contact_area_ && store_contact_data_)
156 {
157 const char* fixarg[10];
158 fixarg[0]="contactAreaConduction";
159 fixarg[1]="all";
160 fixarg[2]="property/atom";
161 fixarg[3]="contactAreaConduction";
162 fixarg[4]="scalar";
163 fixarg[5]="no";
164 fixarg[6]="yes";
165 fixarg[7]="no";
166 fixarg[8]="0.";
167 fix_conduction_contact_area_ = modify->add_fix_property_atom(9,const_cast<char**>(fixarg),style);
168 }
169
170 fix_n_conduction_contacts_ = static_cast<FixPropertyAtom*>(modify->find_fix_property("nContactsConduction","property/atom","scalar",0,0,this->style,false));
171 if(!fix_n_conduction_contacts_ && store_contact_data_)
172 {
173 const char* fixarg[10];
174 fixarg[0]="nContactsConduction";
175 fixarg[1]="all";
176 fixarg[2]="property/atom";
177 fixarg[3]="nContactsConduction";
178 fixarg[4]="scalar";
179 fixarg[5]="no";
180 fixarg[6]="yes";
181 fixarg[7]="no";
182 fixarg[8]="0.";
183 fix_n_conduction_contacts_ = modify->add_fix_property_atom(9,const_cast<char**>(fixarg),style);
184 }
185
186 fix_wall_heattransfer_coeff_ = static_cast<FixPropertyAtom*>(modify->find_fix_property("wallHeattransferCoeff","property/atom","scalar",0,0,this->style,false));
187 if(!fix_wall_heattransfer_coeff_ && store_contact_data_)
188 {
189 const char* fixarg[10];
190 fixarg[0]="wallHeattransferCoeff";
191 fixarg[1]="all";
192 fixarg[2]="property/atom";
193 fixarg[3]="wallHeattransferCoeff";
194 fixarg[4]="scalar";
195 fixarg[5]="no";
196 fixarg[6]="yes";
197 fixarg[7]="no";
198 fixarg[8]="0.";
199 fix_wall_heattransfer_coeff_ = modify->add_fix_property_atom(9,const_cast<char**>(fixarg),style);
200 }
201
202 fix_wall_temperature_ = static_cast<FixPropertyAtom*>(modify->find_fix_property("wallTemp","property/atom","scalar",0,0,this->style,false));
203 if(!fix_wall_temperature_ && store_contact_data_)
204 {
205 const char* fixarg[10];
206 fixarg[0]="wallTemp";
207 fixarg[1]="all";
208 fixarg[2]="property/atom";
209 fixarg[3]="wallTemp";
210 fixarg[4]="scalar";
211 fixarg[5]="no";
212 fixarg[6]="yes";
213 fixarg[7]="no";
214 fixarg[8]="0.";
215 fix_wall_temperature_ = modify->add_fix_property_atom(9,const_cast<char**>(fixarg),style);
216 }
217
218 if(store_contact_data_ && (!fix_conduction_contact_area_ || !fix_n_conduction_contacts_ || !fix_wall_heattransfer_coeff_ || !fix_wall_temperature_))
219 error->one(FLERR,"internal error");
220 }
221
222 /* ---------------------------------------------------------------------- */
223
pre_delete(bool unfixflag)224 void FixHeatGranCond::pre_delete(bool unfixflag)
225 {
226
227 // tell cpl that this fix is deleted
228 if(cpl && unfixflag) cpl->reference_deleted();
229
230 }
231
232 /* ---------------------------------------------------------------------- */
233
setmask()234 int FixHeatGranCond::setmask()
235 {
236 int mask = FixHeatGran::setmask();
237 mask |= PRE_FORCE;
238 mask |= POST_FORCE;
239 return mask;
240 }
241
242 /* ---------------------------------------------------------------------- */
243
updatePtrs()244 void FixHeatGranCond::updatePtrs()
245 {
246 FixHeatGran::updatePtrs();
247
248 if(store_contact_data_)
249 {
250 conduction_contact_area_ = fix_conduction_contact_area_->vector_atom;
251 n_conduction_contacts_ = fix_n_conduction_contacts_->vector_atom;
252 wall_heattransfer_coeff_ = fix_wall_heattransfer_coeff_->vector_atom;
253 wall_temp_ = fix_wall_temperature_->vector_atom;
254 }
255 }
256
257 /* ---------------------------------------------------------------------- */
258
init()259 void FixHeatGranCond::init()
260 {
261 // initialize base class
262 FixHeatGran::init();
263
264 const double *Y, *nu, *Y_orig;
265 double expo, Yeff_ij, Yeff_orig_ij, ratio;
266 int max_type = atom->get_properties()->max_type();
267
268 if (conductivity_) delete []conductivity_;
269 conductivity_ = new double[max_type];
270 fix_conductivity_ =
271 static_cast<FixPropertyGlobal*>(modify->find_fix_property("thermalConductivity","property/global","peratomtype",max_type,0,style));
272
273 // pre-calculate conductivity for possible contact material combinations
274 for(int i=1;i< max_type+1; i++)
275 for(int j=1;j<max_type+1;j++)
276 {
277 conductivity_[i-1] = fix_conductivity_->compute_vector(i-1);
278 if(conductivity_[i-1] < 0.)
279 error->all(FLERR,"Fix heat/gran/conduction: Thermal conductivity must not be < 0");
280 }
281
282 // calculate heat transfer correction
283
284 if(area_correction_flag_)
285 {
286 if(!force->pair_match("gran",0))
287 error->fix_error(FLERR,this,"area correction only works with using granular pair styles");
288
289 expo = 1./pair_gran->stressStrainExponent();
290
291 Y = static_cast<FixPropertyGlobal*>(modify->find_fix_property("youngsModulus","property/global","peratomtype",max_type,0,style))->get_values();
292 nu = static_cast<FixPropertyGlobal*>(modify->find_fix_property("poissonsRatio","property/global","peratomtype",max_type,0,style))->get_values();
293 Y_orig = static_cast<FixPropertyGlobal*>(modify->find_fix_property("youngsModulusOriginal","property/global","peratomtype",max_type,0,style))->get_values();
294
295 // allocate a new array within youngsModulusOriginal
296 static_cast<FixPropertyGlobal*>(modify->find_fix_property("youngsModulusOriginal","property/global","peratomtype",max_type,0,style))->new_array(max_type,max_type);
297
298 // feed deltan_ratio into this array
299 for(int i = 1; i < max_type+1; i++)
300 {
301 for(int j = 1; j < max_type+1; j++)
302 {
303 Yeff_ij = 1./((1.-pow(nu[i-1],2.))/Y[i-1] +(1.-pow(nu[j-1],2.))/Y[j-1]);
304 Yeff_orig_ij = 1./((1.-pow(nu[i-1],2.))/Y_orig[i-1]+(1.-pow(nu[j-1],2.))/Y_orig[j-1]);
305 ratio = pow(Yeff_ij/Yeff_orig_ij,expo);
306
307 static_cast<FixPropertyGlobal*>(modify->find_fix_property("youngsModulusOriginal","property/global","peratomtype",max_type,0,style))->array_modify(i-1,j-1,ratio);
308 }
309 }
310
311 // get reference to deltan_ratio
312 deltan_ratio_ = static_cast<FixPropertyGlobal*>(modify->find_fix_property("youngsModulusOriginal","property/global","peratomtype",max_type,0,style))->get_array_modified();
313 }
314
315 updatePtrs();
316
317 // error checks on coarsegraining
318
319 }
320
321 /* ---------------------------------------------------------------------- */
322
pre_force(int vflag)323 void FixHeatGranCond::pre_force(int vflag)
324 {
325
326 if(store_contact_data_)
327 {
328 fix_wall_heattransfer_coeff_->set_all(0.);
329 fix_wall_temperature_->set_all(0.);
330 }
331 }
332
333 /* ---------------------------------------------------------------------- */
334
post_force(int vflag)335 void FixHeatGranCond::post_force(int vflag)
336 {
337
338 if(history_flag == 0 && CONDUCTION_CONTACT_AREA_OVERLAP == area_calculation_mode_)
339 post_force_eval<0,CONDUCTION_CONTACT_AREA_OVERLAP>(vflag,0);
340 if(history_flag == 1 && CONDUCTION_CONTACT_AREA_OVERLAP == area_calculation_mode_)
341 post_force_eval<1,CONDUCTION_CONTACT_AREA_OVERLAP>(vflag,0);
342
343 if(history_flag == 0 && CONDUCTION_CONTACT_AREA_CONSTANT == area_calculation_mode_)
344 post_force_eval<0,CONDUCTION_CONTACT_AREA_CONSTANT>(vflag,0);
345 if(history_flag == 1 && CONDUCTION_CONTACT_AREA_CONSTANT == area_calculation_mode_)
346 post_force_eval<1,CONDUCTION_CONTACT_AREA_CONSTANT>(vflag,0);
347
348 if(history_flag == 0 && CONDUCTION_CONTACT_AREA_PROJECTION == area_calculation_mode_)
349 post_force_eval<0,CONDUCTION_CONTACT_AREA_PROJECTION>(vflag,0);
350 if(history_flag == 1 && CONDUCTION_CONTACT_AREA_PROJECTION == area_calculation_mode_)
351 post_force_eval<1,CONDUCTION_CONTACT_AREA_PROJECTION>(vflag,0);
352 }
353
354 /* ---------------------------------------------------------------------- */
355
cpl_evaluate(ComputePairGranLocal * caller)356 void FixHeatGranCond::cpl_evaluate(ComputePairGranLocal *caller)
357 {
358 if(caller != cpl) error->all(FLERR,"Illegal situation in FixHeatGranCond::cpl_evaluate");
359
360 if(history_flag == 0 && CONDUCTION_CONTACT_AREA_OVERLAP == area_calculation_mode_)
361 post_force_eval<0,CONDUCTION_CONTACT_AREA_OVERLAP>(0,1);
362 if(history_flag == 1 && CONDUCTION_CONTACT_AREA_OVERLAP == area_calculation_mode_)
363 post_force_eval<1,CONDUCTION_CONTACT_AREA_OVERLAP>(0,1);
364
365 if(history_flag == 0 && CONDUCTION_CONTACT_AREA_CONSTANT == area_calculation_mode_)
366 post_force_eval<0,CONDUCTION_CONTACT_AREA_CONSTANT>(0,1);
367 if(history_flag == 1 && CONDUCTION_CONTACT_AREA_CONSTANT == area_calculation_mode_)
368 post_force_eval<1,CONDUCTION_CONTACT_AREA_CONSTANT>(0,1);
369
370 if(history_flag == 0 && CONDUCTION_CONTACT_AREA_PROJECTION == area_calculation_mode_)
371 post_force_eval<0,CONDUCTION_CONTACT_AREA_PROJECTION>(0,1);
372 if(history_flag == 1 && CONDUCTION_CONTACT_AREA_PROJECTION == area_calculation_mode_)
373 post_force_eval<1,CONDUCTION_CONTACT_AREA_PROJECTION>(0,1);
374 }
375
376 /* ---------------------------------------------------------------------- */
377
378 template <int HISTFLAG,int CONTACTAREA>
post_force_eval(int vflag,int cpl_flag)379 void FixHeatGranCond::post_force_eval(int vflag,int cpl_flag)
380 {
381 double hc,contactArea,delta_n,flux,dirFlux[3];
382 int i,j,ii,jj,inum,jnum;
383 double xtmp,ytmp,ztmp,delx,dely,delz;
384 double radi,radj,radsum,rsq,r,tcoi,tcoj;
385 int *ilist,*jlist,*numneigh,**firstneigh;
386 int *contact_flag,**first_contact_flag;
387
388 int newton_pair = force->newton_pair;
389
390 if (strcmp(force->pair_style,"hybrid")==0)
391 error->warning(FLERR,"Fix heat/gran/conduction implementation may not be valid for pair style hybrid");
392 if (strcmp(force->pair_style,"hybrid/overlay")==0)
393 error->warning(FLERR,"Fix heat/gran/conduction implementation may not be valid for pair style hybrid/overlay");
394
395 inum = pair_gran->list->inum;
396 ilist = pair_gran->list->ilist;
397 numneigh = pair_gran->list->numneigh;
398 firstneigh = pair_gran->list->firstneigh;
399 if(HISTFLAG) first_contact_flag = pair_gran->listgranhistory->firstneigh;
400
401 double *radius = atom->radius;
402 double **x = atom->x;
403 int *type = atom->type;
404 int nlocal = atom->nlocal;
405 int *mask = atom->mask;
406
407 updatePtrs();
408
409 if(store_contact_data_)
410 {
411 fix_conduction_contact_area_->set_all(0.);
412 fix_n_conduction_contacts_->set_all(0.);
413 }
414
415 // loop over neighbors of my atoms
416 for (ii = 0; ii < inum; ii++) {
417 i = ilist[ii];
418 xtmp = x[i][0];
419 ytmp = x[i][1];
420 ztmp = x[i][2];
421 radi = radius[i];
422 jlist = firstneigh[i];
423 jnum = numneigh[i];
424 if(HISTFLAG) contact_flag = first_contact_flag[i];
425
426 for (jj = 0; jj < jnum; jj++) {
427 j = jlist[jj];
428 j &= NEIGHMASK;
429
430 if (!(mask[i] & groupbit) && !(mask[j] & groupbit)) continue;
431
432 if(!HISTFLAG)
433 {
434 delx = xtmp - x[j][0];
435 dely = ytmp - x[j][1];
436 delz = ztmp - x[j][2];
437 rsq = delx*delx + dely*dely + delz*delz;
438 radj = radius[j];
439 radsum = radi + radj;
440 }
441
442 if ((HISTFLAG && contact_flag[jj]) || (!HISTFLAG && (rsq < radsum*radsum))) { //contact
443
444 if(HISTFLAG)
445 {
446 delx = xtmp - x[j][0];
447 dely = ytmp - x[j][1];
448 delz = ztmp - x[j][2];
449 rsq = delx*delx + dely*dely + delz*delz;
450 radj = radius[j];
451 radsum = radi + radj;
452 if(rsq >= radsum*radsum) continue;
453 }
454
455 r = sqrt(rsq);
456
457 if(CONTACTAREA == CONDUCTION_CONTACT_AREA_OVERLAP)
458 {
459
460 if(area_correction_flag_)
461 {
462 delta_n = radsum - r;
463 delta_n *= deltan_ratio_[type[i]-1][type[j]-1];
464 r = radsum - delta_n;
465 }
466
467 if (r < fmax(radi, radj)) // one sphere is inside the other
468 {
469 // set contact area to area of smaller sphere
470 contactArea = fmin(radi,radj);
471 contactArea *= contactArea * M_PI;
472 }
473 else
474 //contact area of the two spheres
475 contactArea = - M_PI/4.0 * ( (r-radi-radj)*(r+radi-radj)*(r-radi+radj)*(r+radi+radj) )/(r*r);
476 }
477 else if (CONTACTAREA == CONDUCTION_CONTACT_AREA_CONSTANT)
478 contactArea = fixed_contact_area_;
479 else if (CONTACTAREA == CONDUCTION_CONTACT_AREA_PROJECTION)
480 {
481 double rmax = std::max(radi,radj);
482 contactArea = M_PI*rmax*rmax;
483 }
484
485 tcoi = conductivity_[type[i]-1];
486 tcoj = conductivity_[type[j]-1];
487 if (tcoi < SMALL_FIX_HEAT_GRAN || tcoj < SMALL_FIX_HEAT_GRAN) hc = 0.;
488 else hc = 4.*tcoi*tcoj/(tcoi+tcoj)*sqrt(contactArea);
489
490 flux = (Temp[j]-Temp[i])*hc;
491
492 dirFlux[0] = flux*delx;
493 dirFlux[1] = flux*dely;
494 dirFlux[2] = flux*delz;
495 if(!cpl_flag)
496 {
497 //Add half of the flux (located at the contact) to each particle in contact
498 heatFlux[i] += flux;
499 directionalHeatFlux[i][0] += 0.50 * dirFlux[0];
500 directionalHeatFlux[i][1] += 0.50 * dirFlux[1];
501 directionalHeatFlux[i][2] += 0.50 * dirFlux[2];
502
503 if(store_contact_data_)
504 {
505 conduction_contact_area_[i] += contactArea;
506 n_conduction_contacts_[i] += 1.;
507 }
508 if (newton_pair || j < nlocal)
509 {
510 heatFlux[j] -= flux;
511 directionalHeatFlux[j][0] += 0.50 * dirFlux[0];
512 directionalHeatFlux[j][1] += 0.50 * dirFlux[1];
513 directionalHeatFlux[j][2] += 0.50 * dirFlux[2];
514
515 if(store_contact_data_)
516 {
517 conduction_contact_area_[j] += contactArea;
518 n_conduction_contacts_[j] += 1.;
519 }
520 }
521 }
522
523 if(cpl_flag && cpl) cpl->add_heat(i,j,flux);
524 }
525 }
526 }
527
528 if(newton_pair)
529 {
530 fix_heatFlux->do_reverse_comm();
531 fix_directionalHeatFlux->do_reverse_comm();
532 fix_conduction_contact_area_->do_reverse_comm();
533 fix_n_conduction_contacts_->do_reverse_comm();
534 }
535
536 if(!cpl_flag && store_contact_data_)
537 for(int i = 0; i < nlocal; i++)
538 {
539 if(n_conduction_contacts_[i] > 0.5)
540 conduction_contact_area_[i] /= n_conduction_contacts_[i];
541 }
542 }
543
544 /* ----------------------------------------------------------------------
545 register and unregister callback to compute
546 ------------------------------------------------------------------------- */
547
register_compute_pair_local(ComputePairGranLocal * ptr)548 void FixHeatGranCond::register_compute_pair_local(ComputePairGranLocal *ptr)
549 {
550
551 if(cpl != NULL)
552 error->all(FLERR,"Fix heat/gran/conduction allows only one compute of type pair/local");
553 cpl = ptr;
554 }
555
unregister_compute_pair_local(ComputePairGranLocal * ptr)556 void FixHeatGranCond::unregister_compute_pair_local(ComputePairGranLocal *ptr)
557 {
558
559 if(cpl != ptr)
560 error->all(FLERR,"Illegal situation in FixHeatGranCond::unregister_compute_pair_local");
561 cpl = NULL;
562 }
563