ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/applications/hydrodynamics/HydrodynamicsModel.cpp
(Generate patch)

Comparing trunk/src/applications/hydrodynamics/HydrodynamicsModel.cpp (file contents):
Revision 893 by tim, Fri Feb 24 21:17:05 2006 UTC vs.
Revision 895 by tim, Mon Mar 13 22:42:40 2006 UTC

# Line 77 | Line 77 | bool HydrodynamicsModel::calcHydrodyanmicsProps() {
77          std::cout << "can not create beads" << std::endl;
78          return false;
79      }
80 +
81 +    //calcResistanceTensor();
82 +    calcDiffusionTensor();
83      
84 + /*    
85      int nbeads = beads_.size();
86      DynamicRectMatrix<double> B(3*nbeads, 3*nbeads);
87      DynamicRectMatrix<double> C(3*nbeads, 3*nbeads);
# Line 134 | Line 138 | bool HydrodynamicsModel::calcHydrodyanmicsProps() {
138  
139      double volume = 0.0;
140      for (std::vector<BeadParam>::iterator iter = beads_.begin(); iter != beads_.end(); ++iter) {
141 <        volume = 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3);
141 >        volume += 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3);
142      }
143          
144      for (std::size_t i = 0; i < nbeads; ++i) {
# Line 144 | Line 148 | bool HydrodynamicsModel::calcHydrodyanmicsProps() {
148              
149              Xitt += Cij;
150              Xitr += U[i] * Cij;
151 <            Xirr += -U[i] * Cij * U[j];            
152 <            //Xirr += -U[i] * Cij * U[j] + (0.166*6 * viscosity_ * volume) * I;            
151 >            //Xirr += -U[i] * Cij * U[j];            
152 >            Xirr += -U[i] * Cij * U[j] + (0.166*6 * viscosity_ * volume) * I;            
153          }
154      }
155  
# Line 230 | Line 234 | bool HydrodynamicsModel::calcHydrodyanmicsProps() {
234      props_.transDiff = Ddtt;
235      props_.transRotDiff = Ddtr;
236      props_.rotDiff = Ddrr;
237 <
237 > */
238      return true;    
239   }
240  
241 < void HydrodynamicsModel::writeBeads(std::ostream& os) {
242 <    std::vector<BeadParam>::iterator iter;
243 <    os << beads_.size() << std::endl;
244 <    os << "Generated by Hydro" << std::endl;
245 <    for (iter = beads_.begin(); iter != beads_.end(); ++iter) {
246 <        os << iter->atomName << "\t" << iter->pos[0] << "\t" << iter->pos[1] << "\t" << iter->pos[2] << std::endl;
241 > void HydrodynamicsModel::calcResistanceTensor() {
242 >    
243 >    int nbeads = beads_.size();
244 >    DynamicRectMatrix<double> B(3*nbeads, 3*nbeads);
245 >    DynamicRectMatrix<double> C(3*nbeads, 3*nbeads);
246 >    Mat3x3d I;
247 >    I(0, 0) = 1.0;
248 >    I(1, 1) = 1.0;
249 >    I(2, 2) = 1.0;
250 >    
251 >    for (std::size_t i = 0; i < nbeads; ++i) {
252 >        for (std::size_t j = 0; j < nbeads; ++j) {
253 >            Mat3x3d Tij;
254 >            if (i != j ) {
255 >                Vector3d Rij = beads_[i].pos - beads_[j].pos;
256 >                double rij = Rij.length();
257 >                double rij2 = rij * rij;
258 >                double sumSigma2OverRij2 = ((beads_[i].radius*beads_[i].radius) + (beads_[i].radius*beads_[i].radius)) / rij2;                
259 >                Mat3x3d tmpMat;
260 >                tmpMat = outProduct(Rij, Rij) / rij2;
261 >                double constant = 8.0 * NumericConstant::PI * viscosity_ * rij;
262 >                Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant;
263 >            }else {
264 >                double constant = 1.0 / (6.0 * NumericConstant::PI * viscosity_ * beads_[i].radius);
265 >                Tij(0, 0) = constant;
266 >                Tij(1, 1) = constant;
267 >                Tij(2, 2) = constant;
268 >            }
269 >            B.setSubMatrix(i*3, j*3, Tij);
270 >        }
271      }
272  
273 < }
273 >  
274 >    //invert B Matrix
275 >    invertMatrix(B, C);
276  
277 < void HydrodynamicsModel::writeDiffCenterAndDiffTensor(std::ostream& os) {
278 <    os << "//viscosity = " << viscosity_ << std::endl;
279 <    os << "//temperature = " << temperature_<< std::endl;
280 <    std::vector<BeadParam>::iterator iter;
281 <    os << sd_->getType() << "\n";
277 >    //prepare U Matrix relative to arbitrary origin O(0.0, 0.0, 0.0)
278 >    std::vector<Mat3x3d> U;
279 >    for (int i = 0; i < nbeads; ++i) {
280 >        Mat3x3d currU;
281 >        currU.setupSkewMat(beads_[i].pos);
282 >        U.push_back(currU);
283 >    }
284 >    
285 >    //calculate Xi matrix at arbitrary origin O
286 >    Mat3x3d Xiott;
287 >    Mat3x3d Xiorr;
288 >    Mat3x3d Xiotr;
289  
290 <    os << "//diffusion center" << std::endl;
254 <    os << props_.diffCenter << std::endl;
290 >    //calculate the total volume
291  
292 <    os << "//translational diffusion tensor" << std::endl;
293 <    os << props_.transDiff << std::endl;
292 >    double volume = 0.0;
293 >    for (std::vector<BeadParam>::iterator iter = beads_.begin(); iter != beads_.end(); ++iter) {
294 >        volume += 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3);
295 >    }
296 >        
297 >    for (std::size_t i = 0; i < nbeads; ++i) {
298 >        for (std::size_t j = 0; j < nbeads; ++j) {
299 >            Mat3x3d Cij;
300 >            C.getSubMatrix(i*3, j*3, Cij);
301 >            
302 >            Xiott += Cij;
303 >            Xiotr += U[i] * Cij;
304 >            //Xiorr += -U[i] * Cij * U[j];            
305 >            Xiorr += -U[i] * Cij * U[j] + (6 * viscosity_ * volume) * I;            
306 >        }
307 >    }
308  
309 <    os << "//translation-rotation coupling diffusion tensor" << std::endl;
310 <    os << props_.transRotDiff << std::endl;
311 <
312 <    os << "//rotational diffusion tensor" << std::endl;
313 <    os << props_.rotDiff << std::endl;
314 <    
315 <    /*
316 <    os << props_.diffCenter[0] << "\t" << props_.diffCenter[1] << "\t" << props_.diffCenter[2] << "\n"
317 <
318 <    os << props_.transDiff(0, 0) << "\t" << props_.transDiff(0, 1) << "\t" << props_.transDiff(0, 2) << "\t"
319 <        << props_.transDiff(1, 0) << "\t" << props_.transDiff(1, 1) << "\t" << props_.transDiff(1, 2) << "\t"
320 <        << props_.transDiff(2, 0) << "\t" << props_.transDiff(2, 1) << "\t" << props_.transDiff(2, 2) << "\n";
321 <    
322 <    os << props_.transRotDiff(0, 0) << "\t" << props_.transRotDiff(0, 1) << "\t" << props_.transRotDiff(0, 2) << "\t"
323 <        << props_.transRotDiff(1, 0) << "\t" << props_.transRotDiff(1, 1) << "\t" << props_.transRotDiff(1, 2) << "\t"
324 <        << props_.transRotDiff(2, 0) << "\t" << props_.transRotDiff(2, 1) << "\t" << props_.transRotDiff(2, 2) << "\t"
309 >    Mat3x3d tmp;
310 >    Mat3x3d tmpInv;
311 >    Vector3d tmpVec;
312 >    tmp(0, 0) = Xiott(1, 1) + Xiott(2, 2);
313 >    tmp(0, 1) = - Xiott(0, 1);
314 >    tmp(0, 2) = -Xiott(0, 2);
315 >    tmp(1, 0) = -Xiott(0, 1);
316 >    tmp(1, 1) = Xiott(0, 0)  + Xiott(2, 2);
317 >    tmp(1, 2) = -Xiott(1, 2);
318 >    tmp(2, 0) = -Xiott(0, 2);
319 >    tmp(2, 1) = -Xiott(1, 2);
320 >    tmp(2, 2) = Xiott(1, 1) + Xiott(0, 0);
321 >    tmpVec[0] = Xiotr(2, 1) - Xiotr(1, 2);
322 >    tmpVec[1] = Xiotr(0, 2) - Xiotr(2, 0);
323 >    tmpVec[2] = Xiotr(1, 0) - Xiotr(0, 1);
324 >    tmpInv = tmp.inverse();    
325 >    Vector3d ror = tmpInv * tmpVec; //center of resistance
326 >    Mat3x3d Uor;
327 >    Uor.setupSkewMat(ror);
328 >    
329 >    Mat3x3d Xirtt;
330 >    Mat3x3d Xirrr;
331 >    Mat3x3d Xirtr;
332  
333 <    os << props_.rotDiff(0, 0) << "\t" << props_.rotDiff(0, 1) << "\t" << props_.rotDiff(0, 2) << "\t"
334 <        << props_.rotDiff(1, 0) << "\t" << props_.rotDiff(1, 1) << "\t" << props_.rotDiff(1, 2) << "\t"
335 <        << props_.rotDiff(2, 0) << "\t" << props_.rotDiff(2, 1) << "\t" << props_.rotDiff(2, 2) << ";"
336 <        << std::endl;
333 >    Xirtt = Xiott;
334 >    Xirtr = (Xiotr - Uor * Xiott) * 1E-8;
335 >    Xirrr = Xiorr - Uor * Xiott * Uor + Xiotr * Uor - Uor * Xiotr.transpose() * 1E-16;
336 > /*
337 >    SquareMatrix<double,6> Xir6x6;
338 >    SquareMatrix<double,6> Dr6x6;
339 >
340 >    Xir6x6.setSubMatrix(0, 0, Xirtt);
341 >    Xir6x6.setSubMatrix(0, 3, Xirtr.transpose());
342 >    Xir6x6.setSubMatrix(3, 0, Xirtr);
343 >    Xir6x6.setSubMatrix(3, 3, Xirrr);
344 >
345 >    invertMatrix(Xir6x6, Dr6x6);
346 >    Mat3x3d Drtt;
347 >    Mat3x3d Drtr;
348 >    Mat3x3d Drrr;
349 >    Dr6x6.getSubMatrix(0, 0, Drtt);
350 >    Dr6x6.getSubMatrix(3, 0, Drtr);
351 >    Dr6x6.getSubMatrix(3, 3, Drrr);
352 >    double kt = OOPSEConstant::kB * temperature_ * 1.66E-2;
353 >    Drtt *= kt;
354 >    Drtr *= kt*1E8;
355 >    Drrr *= kt*1E16;
356 > */  
357 >
358 >    const static Mat3x3d zeroMat(0.0);
359 >
360 >
361 >    
362 >    Mat3x3d XirttInv(0.0);
363 >    XirttInv = Xirtt.inverse();
364 >    
365 >    //Xirr may not be inverted,if it one of the diagonal element is zero, for example
366 >    //( a11 a12 0)
367 >    //( a21 a22 0)
368 >    //( 0    0    0)
369 >    Mat3x3d XirrrInv;
370 >    XirrrInv = Xirrr.inverse();
371 >    tmp = Xirtt - Xirtr.transpose() * XirrrInv * Xirtr;
372 >    tmpInv = tmp.inverse();
373 >
374 >    Mat3x3d Drtt;
375 >    Mat3x3d Drtr;
376 >    Mat3x3d Drrr;
377 >    double kt = OOPSEConstant::kB * temperature_ * 1.66E-2;
378 >    Drtt = kt * tmpInv;
379 >    Drtr = -kt*XirrrInv * Xirtr * tmpInv* 1.0E8;
380 >
381 >    tmp = Xirrr - Xirtr * XirttInv * Xirtr.transpose();    
382 >    tmpInv = tmp.inverse();
383 >    
384 >    Drrr = kt * tmpInv*1.0E16;
385 >
386 >    std::cout << "-----------------------------------------\n";
387 >    std::cout << "center of resistance :" << std::endl;
388 >    std::cout << ror << std::endl;
389 >    std::cout << "resistant tensor at center of resistance" << std::endl;
390 >    std::cout << "translation:" << std::endl;
391 >    std::cout << Xirtt << std::endl;
392 >    std::cout << "translation-rotation:" << std::endl;
393 >    std::cout << Xirtr << std::endl;
394 >    std::cout << "rotation:" << std::endl;
395 >    std::cout << Xirrr << std::endl;
396 >    std::cout << "diffusion tensor at center of resistance" << std::endl;
397 >    std::cout << "translation:" << std::endl;
398 >    std::cout << Drtt << std::endl;
399 >    std::cout << "translation-rotation:" << std::endl;
400 >    std::cout << Drtr << std::endl;
401 >    std::cout << "rotation:" << std::endl;
402 >    std::cout << Drrr << std::endl;    
403 >    std::cout << "-----------------------------------------\n";
404 >
405 > }
406 >
407 > void HydrodynamicsModel::calcDiffusionTensor() {
408 >    int nbeads = beads_.size();
409 >    DynamicRectMatrix<double> B(3*nbeads, 3*nbeads);
410 >    DynamicRectMatrix<double> C(3*nbeads, 3*nbeads);
411 >    Mat3x3d I;
412 >    I(0, 0) = 1.0;
413 >    I(1, 1) = 1.0;
414 >    I(2, 2) = 1.0;
415 >    
416 >    for (std::size_t i = 0; i < nbeads; ++i) {
417 >        for (std::size_t j = 0; j < nbeads; ++j) {
418 >            Mat3x3d Tij;
419 >            if (i != j ) {
420 >                Vector3d Rij = beads_[i].pos - beads_[j].pos;
421 >                double rij = Rij.length();
422 >                double rij2 = rij * rij;
423 >                double sumSigma2OverRij2 = ((beads_[i].radius*beads_[i].radius) + (beads_[i].radius*beads_[i].radius)) / rij2;                
424 >                Mat3x3d tmpMat;
425 >                tmpMat = outProduct(Rij, Rij) / rij2;
426 >                double constant = 8.0 * NumericConstant::PI * viscosity_ * rij;
427 >                Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant;
428 >            }else {
429 >                double constant = 1.0 / (6.0 * NumericConstant::PI * viscosity_ * beads_[i].radius);
430 >                Tij(0, 0) = constant;
431 >                Tij(1, 1) = constant;
432 >                Tij(2, 2) = constant;
433 >            }
434 >            B.setSubMatrix(i*3, j*3, Tij);
435 >        }
436 >    }
437 >
438 >    //invert B Matrix
439 >    invertMatrix(B, C);
440 >
441 >    //prepare U Matrix relative to arbitrary origin O(0.0, 0.0, 0.0)
442 >    std::vector<Mat3x3d> U;
443 >    for (int i = 0; i < nbeads; ++i) {
444 >        Mat3x3d currU;
445 >        currU.setupSkewMat(beads_[i].pos);
446 >        U.push_back(currU);
447 >    }
448 >    
449 >    //calculate Xi matrix at arbitrary origin O
450 >    Mat3x3d Xitt;
451 >    Mat3x3d Xirr;
452 >    Mat3x3d Xitr;
453 >
454 >    //calculate the total volume
455 >
456 >    double volume = 0.0;
457 >    for (std::vector<BeadParam>::iterator iter = beads_.begin(); iter != beads_.end(); ++iter) {
458 >        volume += 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3);
459 >    }
460 >        
461 >    for (std::size_t i = 0; i < nbeads; ++i) {
462 >        for (std::size_t j = 0; j < nbeads; ++j) {
463 >            Mat3x3d Cij;
464 >            C.getSubMatrix(i*3, j*3, Cij);
465 >            
466 >            Xitt += Cij;
467 >            Xitr += U[i] * Cij;
468 >            //Xirr += -U[i] * Cij * U[j];            
469 >            Xirr += -U[i] * Cij * U[j] + (6 * viscosity_ * volume) * I;            
470 >        }
471 >    }
472 >
473 >    //invert Xi to get Diffusion Tensor at arbitrary origin O
474 >    RectMatrix<double, 6, 6> Xi;    
475 >    RectMatrix<double, 6, 6> Do;
476 >    Xi.setSubMatrix(0, 0, Xitt);
477 >    Xi.setSubMatrix(0, 3, Xitr.transpose());
478 >    Xi.setSubMatrix(3, 0, Xitr);
479 >    Xi.setSubMatrix(3, 3, Xirr);
480 >    //invertMatrix(Xi, Do);
481 >    //double kt = OOPSEConstant::kB * temperature_ * 1.66E-2;
482 >
483 >    //1 poise = 0.1 N.S/m^2 = 1.661E-3 amu/ (Angstrom*fs)
484 >    double kt = OOPSEConstant::kB * temperature_ * 1.66E-3;
485 >
486 >    Mat3x3d Dott; //translational diffusion tensor at arbitrary origin O
487 >    Mat3x3d Dorr; //rotational diffusion tensor at arbitrary origin O
488 >    Mat3x3d Dotr; //translation-rotation couplingl diffusion tensor at arbitrary origin O
489 >
490 >    const static Mat3x3d zeroMat(0.0);
491 >    
492 >    Mat3x3d XittInv(0.0);
493 >    XittInv = Xitt.inverse();
494 >    
495 >    //Xirr may not be inverted,if it one of the diagonal element is zero, for example
496 >    //( a11 a12 0)
497 >    //( a21 a22 0)
498 >    //( 0    0    0)
499 >    Mat3x3d XirrInv;
500 >    XirrInv = Xirr.inverse();
501 >
502 >    Mat3x3d tmp;
503 >    Mat3x3d tmpInv;
504 >    tmp = Xitt - Xitr.transpose() * XirrInv * Xitr;
505 >    tmpInv = tmp.inverse();
506 >
507 >    //Dott = kt * tmpInv; //unit in A^2/fs
508 >    Dott = tmpInv;
509 >    //Dotr = -kt*XirrInv * Xitr * tmpInv*1E8;
510 >    //Dotr = -kt*XirrInv * Xitr * tmpInv;
511 >    Dotr = -XirrInv* Xitr * tmpInv;
512 >    
513 >    tmp = Xirr - Xitr * XittInv * Xitr.transpose();    
514 >    tmpInv = tmp.inverse();
515 >    
516 >    //Dorr = kt * tmpInv*1E16;
517 >    //Dorr = kt * tmpInv;
518 >    Dorr = tmpInv;
519 >    //calculate center of diffusion
520 >    tmp(0, 0) = Dorr(1, 1) + Dorr(2, 2);
521 >    tmp(0, 1) = - Dorr(0, 1);
522 >    tmp(0, 2) = -Dorr(0, 2);
523 >    tmp(1, 0) = -Dorr(0, 1);
524 >    tmp(1, 1) = Dorr(0, 0)  + Dorr(2, 2);
525 >    tmp(1, 2) = -Dorr(1, 2);
526 >    tmp(2, 0) = -Dorr(0, 2);
527 >    tmp(2, 1) = -Dorr(1, 2);
528 >    tmp(2, 2) = Dorr(1, 1) + Dorr(0, 0);
529 >
530 >    Vector3d tmpVec;
531 >    tmpVec[0] = Dotr(1, 2) - Dotr(2, 1);
532 >    tmpVec[1] = Dotr(2, 0) - Dotr(0, 2);
533 >    tmpVec[2] = Dotr(0, 1) - Dotr(1, 0);
534 >
535 >    tmpInv = tmp.inverse();
536 >    
537 >    Vector3d rod = tmpInv * tmpVec;
538 >
539 >    //calculate Diffusion Tensor at center of diffusion
540 >    Mat3x3d Uod;
541 >    Uod.setupSkewMat(rod);
542 >    
543 >    Mat3x3d Ddtt; //translational diffusion tensor at diffusion center
544 >    Mat3x3d Ddtr; //rotational diffusion tensor at diffusion center
545 >    Mat3x3d Ddrr; //translation-rotation couplingl diffusion tensor at diffusion tensor
546 >    
547 >    Ddtt = Dott - Uod * Dorr * Uod + Dotr.transpose() * Uod - Uod * Dotr;
548 >    Ddrr = Dorr;
549 >    Ddtr = Dotr + Dorr * Uod;
550 >
551 >    props_.diffCenter = rod;
552 >    props_.Ddtt = Ddtt;
553 >    props_.Ddtr = Ddtr;
554 >    props_.Ddrr = Ddrr;
555 >
556 >    SquareMatrix<double, 6> Dd;
557 >    Dd.setSubMatrix(0, 0, Ddtt);
558 >    Dd.setSubMatrix(0, 3, Ddtr.transpose());
559 >    Dd.setSubMatrix(3, 0, Ddtr);
560 >    Dd.setSubMatrix(3, 3, Ddrr);    
561 >    SquareMatrix<double, 6> Xid;
562 >    invertMatrix(Dd, Xid);
563 >
564 >    Ddtt *= kt;
565 >    Ddtr *= kt;
566 >    Ddrr *= kt;
567 >    Xid /= 1.66E-3;
568 >
569 >    Xid.getSubMatrix(0, 0, props_.Xidtt);
570 >    Xid.getSubMatrix(0, 3, props_.Xidrt);
571 >    Xid.getSubMatrix(3, 0, props_.Xidtr);
572 >    Xid.getSubMatrix(3, 3, props_.Xidrr);
573 >
574 >    /*        
575 >    std::cout << "center of diffusion :" << std::endl;
576 >    std::cout << rod << std::endl;
577 >    std::cout << "diffusion tensor at center of diffusion" << std::endl;
578 >    std::cout << "translation:" << std::endl;
579 >    std::cout << Ddtt << std::endl;
580 >    std::cout << "translation-rotation:" << std::endl;
581 >    std::cout << Ddtr << std::endl;
582 >    std::cout << "rotation:" << std::endl;
583 >    std::cout << Ddrr << std::endl;
584      */
585 +      
586   }
587  
588 + void HydrodynamicsModel::writeBeads(std::ostream& os) {
589 +    std::vector<BeadParam>::iterator iter;
590 +    os << beads_.size() << std::endl;
591 +    os << "Generated by Hydro" << std::endl;
592 +    for (iter = beads_.begin(); iter != beads_.end(); ++iter) {
593 +        os << iter->atomName << "\t" << iter->pos[0] << "\t" << iter->pos[1] << "\t" << iter->pos[2] << std::endl;
594 +    }
595 +
596   }
597 +
598 + void HydrodynamicsModel::writeDiffCenterAndDiffTensor(std::ostream& os) {
599 +
600 +    os << sd_->getType() << "\t";
601 +    os << props_.diffCenter[0] << "\t" << props_.diffCenter[1] << "\t" << props_.diffCenter[2] << "\t";
602 +
603 +    os << props_.Ddtt(0, 0) << "\t" << props_.Ddtt(0, 1) << "\t" << props_.Ddtt(0, 2) << "\t"
604 +        << props_.Ddtt(1, 0) << "\t" << props_.Ddtt(1, 1) << "\t" << props_.Ddtt(1, 2) << "\t"
605 +        << props_.Ddtt(2, 0) << "\t" << props_.Ddtt(2, 1) << "\t" << props_.Ddtt(2, 2) << "\t";
606 +    
607 +    os << props_.Ddtr(0, 0) << "\t" << props_.Ddtr(0, 1) << "\t" << props_.Ddtr(0, 2) << "\t"
608 +        << props_.Ddtr(1, 0) << "\t" << props_.Ddtr(1, 1) << "\t" << props_.Ddtr(1, 2) << "\t"
609 +        << props_.Ddtr(2, 0) << "\t" << props_.Ddtr(2, 1) << "\t" << props_.Ddtr(2, 2) << "\t";
610 +
611 +    os << props_.Ddrr(0, 0) << "\t" << props_.Ddrr(0, 1) << "\t" << props_.Ddrr(0, 2) << "\t"
612 +        << props_.Ddrr(1, 0) << "\t" << props_.Ddrr(1, 1) << "\t" << props_.Ddrr(1, 2) << "\t"
613 +        << props_.Ddrr(2, 0) << "\t" << props_.Ddrr(2, 1) << "\t" << props_.Ddrr(2, 2) <<"\t";        
614 +
615 +    os << props_.Xidtt(0, 0) << "\t" << props_.Xidtt(0, 1) << "\t" << props_.Xidtt(0, 2) << "\t"
616 +        << props_.Xidtt(1, 0) << "\t" << props_.Xidtt(1, 1) << "\t" << props_.Xidtt(1, 2) << "\t"
617 +        << props_.Xidtt(2, 0) << "\t" << props_.Xidtt(2, 1) << "\t" << props_.Xidtt(2, 2) << "\t";
618 +
619 +    os << props_.Xidrt(0, 0) << "\t" << props_.Xidrt(0, 1) << "\t" << props_.Xidrt(0, 2) << "\t"
620 +        << props_.Xidrt(1, 0) << "\t" << props_.Xidrt(1, 1) << "\t" << props_.Xidrt(1, 2) << "\t"
621 +        << props_.Xidrt(2, 0) << "\t" << props_.Xidrt(2, 1) << "\t" << props_.Xidrt(2, 2) << "\t";
622 +    
623 +    os << props_.Xidtr(0, 0) << "\t" << props_.Xidtr(0, 1) << "\t" << props_.Xidtr(0, 2) << "\t"
624 +        << props_.Xidtr(1, 0) << "\t" << props_.Xidtr(1, 1) << "\t" << props_.Xidtr(1, 2) << "\t"
625 +        << props_.Xidtr(2, 0) << "\t" << props_.Xidtr(2, 1) << "\t" << props_.Xidtr(2, 2) << "\t";
626 +
627 +    os << props_.Xidrr(0, 0) << "\t" << props_.Xidrr(0, 1) << "\t" << props_.Xidrr(0, 2) << "\t"
628 +        << props_.Xidrr(1, 0) << "\t" << props_.Xidrr(1, 1) << "\t" << props_.Xidrr(1, 2) << "\t"
629 +        << props_.Xidrr(2, 0) << "\t" << props_.Xidrr(2, 1) << "\t" << props_.Xidrr(2, 2) << std::endl;
630 +    
631 + }
632 +
633 + }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines