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