OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
gnuplotexportmodule.C
Go to the documentation of this file.
1 /*
2  *
3  * ##### ##### ###### ###### ### ###
4  * ## ## ## ## ## ## ## ### ##
5  * ## ## ## ## #### #### ## # ##
6  * ## ## ## ## ## ## ## ##
7  * ## ## ## ## ## ## ## ##
8  * ##### ##### ## ###### ## ##
9  *
10  *
11  * OOFEM : Object Oriented Finite Element Code
12  *
13  * Copyright (C) 1993 - 2013 Borek Patzak
14  *
15  *
16  *
17  * Czech Technical University, Faculty of Civil Engineering,
18  * Department of Structural Mechanics, 166 29 Prague, Czech Republic
19  *
20  * This library is free software; you can redistribute it and/or
21  * modify it under the terms of the GNU Lesser General Public
22  * License as published by the Free Software Foundation; either
23  * version 2.1 of the License, or (at your option) any later version.
24  *
25  * This program is distributed in the hope that it will be useful,
26  * but WITHOUT ANY WARRANTY; without even the implied warranty of
27  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
28  * Lesser General Public License for more details.
29  *
30  * You should have received a copy of the GNU Lesser General Public
31  * License along with this library; if not, write to the Free Software
32  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
33  */
34 
35 /*
36  * gnuplotexportmodule.C
37  *
38  * Created on: Jan 29, 2014
39  * Author: svennine
40  */
41 
42 #include "gnuplotexportmodule.h"
43 #include "classfactory.h"
44 #include "valuemodetype.h"
45 #include "../sm/EngineeringModels/structengngmodel.h"
46 #include "outputmanager.h"
47 #include "dofmanager.h"
48 #include "boundarycondition.h"
49 #include "xfem/enrichmentitem.h"
50 #include "xfem/xfemmanager.h"
51 #include "../sm/Materials/InterfaceMaterials/structuralinterfacematerialstatus.h"
52 #include "xfem/XFEMDebugTools.h"
53 #include "prescribedgradient.h"
56 #include "gausspoint.h"
57 #include "timestep.h"
58 #include "feinterpol.h"
60 #include "dofmanager.h"
61 #include "xfem/matforceevaluator.h"
62 #include "function.h"
65 
66 #include <sstream>
67 
68 namespace oofem {
70 
72  ExportModule(n, e),
73  mExportReactionForces(false),
74  mExportBoundaryConditions(false),
75  mExportBoundaryConditionsExtra(false),
76  mExportMesh(false),
77  mExportXFEM(false),
78  mExportCrackLength(false),
79  mExportInterfaceEl(false),
80  mMonitorNodeIndex(-1),
81  mpMatForceEvaluator( new MaterialForceEvaluator() )
82 {}
83 
85 {}
86 
88 {
96 
98 
100 // printf("mMatForceRadii: "); mMatForceRadii.printYourself();
101 
102  return ExportModule::initializeFrom(ir);
103 }
104 
105 void GnuplotExportModule::doOutput(TimeStep *tStep, bool forcedOutput)
106 {
107  if (!(testTimeStepOutput(tStep) || forcedOutput)) {
108  return;
109  }
110 
111  // Export the sum of reaction forces for each Dirichlet BC
113  outputReactionForces(tStep);
114  }
115 
116  Domain *domain = emodel->giveDomain(1);
117 
118  // Export output from boundary conditions
120  int numBC = domain->giveNumberOfBoundaryConditions();
121 
122  for(int i = 1; i <= numBC; i++) {
123 
124  PrescribedGradient *presGradBC = dynamic_cast<PrescribedGradient*>( domain->giveBc(i) );
125  if(presGradBC != NULL) {
126  outputBoundaryCondition(*presGradBC, tStep);
127  }
128 
129 
130  PrescribedGradientBCNeumann *presGradBCNeumann = dynamic_cast<PrescribedGradientBCNeumann*>( domain->giveBc(i) );
131  if(presGradBCNeumann != NULL) {
132  outputBoundaryCondition(*presGradBCNeumann, tStep);
133  }
134 
135  PrescribedGradientBCWeak *presGradBCWeak = dynamic_cast<PrescribedGradientBCWeak*>( domain->giveBc(i) );
136  if(presGradBCWeak != NULL) {
137  outputBoundaryCondition(*presGradBCWeak, tStep);
138  }
139 
140  }
141  }
142 
143  mTimeHist.push_back( tStep->giveTargetTime() );
144 
145  if(mExportXFEM) {
146  if(domain->hasXfemManager()) {
147  XfemManager *xMan = domain->giveXfemManager();
148 
149  int numEI = xMan->giveNumberOfEnrichmentItems();
150 
151  std::vector< std::vector<FloatArray> > points;
152 
153  for(int i = 1; i <= numEI; i++) {
154  EnrichmentItem *ei = xMan->giveEnrichmentItem(i);
155  ei->callGnuplotExportModule(*this, tStep);
156 
157  GeometryBasedEI *geoEI = dynamic_cast<GeometryBasedEI*>(ei);
158  if(geoEI != NULL) {
159  std::vector<FloatArray> eiPoints;
160  geoEI->giveSubPolygon(eiPoints, 0.0, 1.0);
161  points.push_back(eiPoints);
162  }
163  }
164 
165  outputXFEMGeometry(points);
166  }
167  }
168 
169  if(mExportMesh) {
170  outputMesh(*domain);
171  }
172 
173  if(mMonitorNodeIndex != -1) {
175  outputNodeDisp(*dMan, tStep);
176  }
177 
178  if(mExportInterfaceEl) {
179  outputInterfaceEl(*domain, tStep);
180  }
181 }
182 
184 {
186 }
187 
189 {
190 
191 }
192 
194 // Help functions
196 {
197  // Add sum of reaction forces to arrays
198  // Compute sum of reaction forces for each BC number
199  Domain *domain = emodel->giveDomain(1);
200  StructuralEngngModel *seMod = dynamic_cast<StructuralEngngModel* >(emodel);
201  if(seMod == NULL) {
202  OOFEM_ERROR("failed to cast to StructuralEngngModel.");
203  }
204 
205  FloatArray reactions;
206  IntArray dofManMap, dofidMap, eqnMap;
207 
208  // test if solution step output is active
209  if ( !testTimeStepOutput(tStep) ) {
210  return;
211  }
212 
213  // map contains corresponding dofmanager and dofs numbers corresponding to prescribed equations
214  // sorted according to dofmanger number and as a minor crit. according to dof number
215  // this is necessary for extractor, since the sorted output is expected
216  seMod->buildReactionTable(dofManMap, dofidMap, eqnMap, tStep, 1);
217 
218  // compute reaction forces
219  seMod->computeReaction(reactions, tStep, 1);
220 
221  // Find highest index of prescribed dofs
222  int maxIndPresDof = 0;
223  for ( int i = 1; i <= dofManMap.giveSize(); i++ ) {
224  maxIndPresDof = std::max(maxIndPresDof, dofidMap.at(i));
225  }
226 
227  int numBC = domain->giveNumberOfBoundaryConditions();
228 
229  while ( mReactionForceHistory.size() < size_t(numBC) ) {
230  std::vector<FloatArray> emptyArray;
231  mReactionForceHistory.push_back( emptyArray );
232  }
233 
234  maxIndPresDof = domain->giveNumberOfSpatialDimensions();
235 
236  while ( mDispHist.size() < size_t(numBC) ) {
237  std::vector<double> emptyArray;
238  mDispHist.push_back( emptyArray );
239  }
240 
241  for(int bcInd = 0; bcInd < numBC; bcInd++) {
242  FloatArray fR(maxIndPresDof), disp(numBC);
243  fR.zero();
244 
245 
246  for ( int i = 1; i <= dofManMap.giveSize(); i++ ) {
247  DofManager *dMan = domain->giveDofManager( dofManMap.at(i) );
248  Dof *dof = dMan->giveDofWithID( dofidMap.at(i) );
249 
250  if ( dof->giveBcId() == bcInd+1 ) {
251  fR.at( dofidMap.at(i) ) += reactions.at( eqnMap.at(i) );
252 
253  // Slightly dirty
254  BoundaryCondition *bc = dynamic_cast<BoundaryCondition*> (domain->giveBc(bcInd+1));
255  if ( bc != NULL ) {
256  disp.at(bcInd+1) = std::max( disp.at(bcInd+1), bc->give(dof, VM_Total, tStep->giveTargetTime()) );
257  }
259  }
260  }
261 
262  mDispHist[bcInd].push_back(disp.at(bcInd+1));
263  mReactionForceHistory[bcInd].push_back(fR);
264 
265 
266 
267  // X
268  FILE * pFileX;
269  char fileNameX[100];
270  sprintf(fileNameX, "ReactionForceGnuplotBC%dX.dat", bcInd+1);
271  pFileX = fopen ( fileNameX , "wb" );
272 
273  fprintf(pFileX, "#u Fx\n");
274  for ( size_t j = 0; j < mDispHist[bcInd].size(); j++ ) {
275  fprintf(pFileX, "%e %e\n", mDispHist[bcInd][j], mReactionForceHistory[bcInd][j].at(1) );
276  }
277 
278  fclose(pFileX);
279 
280  // Y
281  FILE * pFileY;
282  char fileNameY[100];
283  sprintf(fileNameY, "ReactionForceGnuplotBC%dY.dat", bcInd+1);
284  pFileY = fopen ( fileNameY , "wb" );
285 
286  fprintf(pFileY, "#u Fx\n");
287  for ( size_t j = 0; j < mDispHist[bcInd].size(); j++ ) {
288  if( mReactionForceHistory[bcInd][j].giveSize() >= 2 ) {
289  fprintf(pFileY, "%e %e\n", mDispHist[bcInd][j], mReactionForceHistory[bcInd][j].at(2) );
290  }
291  }
292 
293  fclose(pFileY);
294 
295  }
296 }
297 
299 {
300 
301 }
302 
304 {
305  const std::vector<GaussPoint*> &czGaussPoints = iCrack.giveCohesiveZoneGaussPoints();
306 
307  std::vector<double> arcLengthPositions, normalJumps, tangJumps, normalTractions, tangTractions;
308 
309  std::vector<double> arcLengthPositionsB = iCrack.giveCohesiveZoneArcPositions();
310 
311  const BasicGeometry *bg = iCrack.giveGeometry();
312 
313 // for( auto *gp: czGaussPoints ) {
314 //
315 // StructuralInterfaceMaterialStatus *matStat = dynamic_cast<StructuralInterfaceMaterialStatus*> ( gp->giveMaterialStatus() );
316 // if(matStat != NULL) {
317 //
318 // // Compute arc length position of the Gauss point
319 // const FloatArray &coord = (gp->giveGlobalCoordinates());
320 // printf("coord: "); coord.printYourself();
321 // double tangDist = 0.0, arcPos = 0.0;
322 // bg->computeTangentialSignDist(tangDist, coord, arcPos);
323 // arcLengthPositions.push_back(arcPos);
324 //
338 // }
339 // }
340 
341 #if 1
342  size_t num_cz_gp = czGaussPoints.size();
343 
344 // for( auto *gp: czGaussPoints ) {
345  for(size_t gp_ind = 0; gp_ind < num_cz_gp; gp_ind++) {
346 // printf("gp_ind: %lu\n", gp_ind );
347  GaussPoint *gp = czGaussPoints[gp_ind];
348 
349  if( gp != NULL ) {
350 
352  if(matStat != NULL) {
353 
354  // Compute arc length position of the Gauss point
355  const FloatArray &coord = (gp->giveGlobalCoordinates());
356  double tangDist = 0.0, arcPos = 0.0;
357  bg->computeTangentialSignDist(tangDist, coord, arcPos);
358 // printf("arcPos: %e\n", arcPos );
359  arcLengthPositions.push_back(arcPos);
360 
361  // Compute displacement jump in normal and tangential direction
362  // Local numbering: (tang_z, tang, normal)
363  const FloatArray &jumpLoc = matStat->giveJump();
364 
365  double normalJump = jumpLoc.at(3);
366  normalJumps.push_back(normalJump);
367 
368 
369  tangJumps.push_back( jumpLoc.at(2) );
370 
371 
372  const FloatArray &trac = matStat->giveFirstPKTraction();
373  normalTractions.push_back(trac.at(3));
374 
375  tangTractions.push_back(trac.at(2));
376  }
377  else {
379 
380  if(fe2ms != NULL) {
381 // printf("Casted to StructuralFE2MaterialStatus.\n");
382 
383  const FloatArray &coord = (gp->giveGlobalCoordinates());
384  double tangDist = 0.0, arcPos = 0.0;
385  bg->computeTangentialSignDist(tangDist, coord, arcPos);
386  // printf("arcPos: %e\n", arcPos );
387  arcLengthPositions.push_back(arcPos);
388 
389  const FloatArray &n = fe2ms->giveNormal();
390 // printf("n: "); n.printYourself();
391 
392  const FloatArray &sig_v = fe2ms->giveStressVector();
393 // printf("sig_v: "); sig_v.printYourself();
394 
395  FloatMatrix sig_m(2,2);
396  sig_m(0,0) = sig_v(0);
397  sig_m(1,1) = sig_v(1);
398  sig_m(0,1) = sig_m(1,0) = sig_v( sig_v.giveSize()-1 );
399 
400  FloatArray trac(2);
401  trac(0) = sig_m(0,0)*n(0) + sig_m(0,1)*n(1);
402  trac(1) = sig_m(1,0)*n(0) + sig_m(1,1)*n(1);
403  double trac_n = trac(0)*n(0) + trac(1)*n(1);
404  normalTractions.push_back(trac_n);
405 
406  const FloatArray t = {-n(1), n(0)};
407  double trac_t = trac(0)*t(0) + trac(1)*t(1);
408  tangTractions.push_back(trac_t);
409 
410 
411  }
412 
413  }
414  }
415  }
416 #endif
417 
418 
419 
420  Domain *domain = emodel->giveDomain(1);
421  XfemManager *xMan = domain->giveXfemManager();
422  if ( xMan != NULL ) {
423  double time = 0.0;
424 
426  if ( ts != NULL ) {
427  time = ts->giveTargetTime();
428  }
429 
430  int eiIndex = iCrack.giveNumber();
431 
432  std :: stringstream strNormalJump;
433  strNormalJump << "NormalJumpGnuplotEI" << eiIndex << "Time" << time << ".dat";
434  std :: string nameNormalJump = strNormalJump.str();
435  XFEMDebugTools::WriteArrayToGnuplot(nameNormalJump, arcLengthPositions, normalJumps);
436 
437  std :: stringstream strTangJump;
438  strTangJump << "TangJumpGnuplotEI" << eiIndex << "Time" << time << ".dat";
439  std :: string nameTangJump = strTangJump.str();
440  XFEMDebugTools::WriteArrayToGnuplot(nameTangJump, arcLengthPositions, tangJumps);
441 
442  std :: stringstream strNormalTrac;
443  strNormalTrac << "NormalTracGnuplotEI" << eiIndex << "Time" << time << ".dat";
444  std :: string nameNormalTrac = strNormalTrac.str();
445  XFEMDebugTools::WriteArrayToGnuplot(nameNormalTrac, arcLengthPositions, normalTractions);
446 
447  std :: stringstream strTangTrac;
448  strTangTrac << "TangTracGnuplotEI" << eiIndex << "Time" << time << ".dat";
449  std :: string nameTangTrac = strTangTrac.str();
450  XFEMDebugTools::WriteArrayToGnuplot(nameTangTrac, arcLengthPositions, tangTractions);
451 
452 
453  std::vector<FloatArray> matForcesStart, matForcesEnd;
454  std::vector<double> radii;
455 
456  // Material forces
457  for(double matForceRadius : mMatForceRadii) {
458 
459  radii.push_back(matForceRadius);
460 
461  EnrichmentFront *efStart = iCrack.giveEnrichmentFrontStart();
462  const TipInfo &tipInfoStart = efStart->giveTipInfo();
463 
464  FloatArray matForceStart;
465  mpMatForceEvaluator->computeMaterialForce(matForceStart, *domain, tipInfoStart, tStep, matForceRadius);
466 
467  if(matForceStart.giveSize() > 0) {
468  matForcesStart.push_back(matForceStart);
469  }
470  else {
471  matForcesStart.push_back({0.0,0.0});
472  }
473 
474 
475  EnrichmentFront *efEnd = iCrack.giveEnrichmentFrontEnd();
476  const TipInfo &tipInfoEnd = efEnd->giveTipInfo();
477 
478  FloatArray matForceEnd;
479  mpMatForceEvaluator->computeMaterialForce(matForceEnd, *domain, tipInfoEnd, tStep, matForceRadius);
480 
481  if(matForceEnd.giveSize() > 0) {
482  matForcesEnd.push_back(matForceEnd);
483  }
484  else {
485  matForcesEnd.push_back({0.0,0.0});
486  }
487 
488  }
489 
490  std::vector< std::vector<FloatArray> > matForcesStartArray, matForcesEndArray;
491  matForcesStartArray.push_back(matForcesStart);
492  matForcesEndArray.push_back(matForcesEnd);
493 
494 
495  std :: stringstream strRadii;
496  strRadii << "MatForceRadiiGnuplotTime" << time << "Crack" << iCrack.giveNumber() << ".dat";
497  XFEMDebugTools::WriteArrayToGnuplot(strRadii.str(), radii, radii);
498 
499 
500  std :: stringstream strMatForcesStart;
501  strMatForcesStart << "MatForcesStartGnuplotTime" << time << "Crack" << iCrack.giveNumber() << ".dat";
502  WritePointsToGnuplot(strMatForcesStart.str(), matForcesStartArray);
503 
504 
505  std :: stringstream strMatForcesEnd;
506  strMatForcesEnd << "MatForcesEndGnuplotTime" << time << "Crack" << iCrack.giveNumber() << ".dat";
507  WritePointsToGnuplot(strMatForcesEnd.str(), matForcesEndArray);
508 
509  double crackLength = iCrack.computeLength();
510  // printf("crackLength: %e\n", crackLength );
511  mCrackLengthHist[eiIndex].push_back(crackLength);
512  std :: stringstream strCrackLength;
513  strCrackLength << "CrackLengthGnuplotEI" << eiIndex << "Time" << time << ".dat";
514  XFEMDebugTools::WriteArrayToGnuplot(strCrackLength.str(), mTimeHist, mCrackLengthHist[eiIndex]);
515  }
516 }
517 
518 void GnuplotExportModule::outputXFEMGeometry(const std::vector< std::vector<FloatArray> > &iEnrItemPoints)
519 {
520  double time = 0.0;
521 
523  if ( ts != NULL ) {
524  time = ts->giveTargetTime();
525  }
526 
527  std :: stringstream strCracks;
528  strCracks << "CracksTime" << time << ".dat";
529  std :: string nameCracks = strCracks.str();
530  WritePointsToGnuplot(nameCracks, iEnrItemPoints);
531 }
532 
534 {
535  FloatArray stress;
536  iBC.computeField(stress, tStep);
537  printf("Mean stress computed in Gnuplot export module: "); stress.printYourself();
538 
539  double time = 0.0;
540 
542  if ( ts != NULL ) {
543  time = ts->giveTargetTime();
544  }
545 
546  int bcIndex = iBC.giveNumber();
547 
548  // Homogenized stress
549  std :: stringstream strMeanStress;
550  strMeanStress << "PrescribedGradientGnuplotMeanStress" << bcIndex << "Time" << time << ".dat";
551  std :: string nameMeanStress = strMeanStress.str();
552  std::vector<double> componentArray, stressArray;
553 
554  for(int i = 1; i <= stress.giveSize(); i++) {
555  componentArray.push_back(i);
556  stressArray.push_back(stress.at(i));
557  }
558 
559  XFEMDebugTools::WriteArrayToGnuplot(nameMeanStress, componentArray, stressArray);
560 
561  FloatArray grad;
562  iBC.giveGradientVoigt(grad);
563  outputGradient(iBC.giveNumber(), *iBC.giveDomain(), grad, tStep);
564 }
565 
567 {
568  FloatArray stress;
569  iBC.computeField(stress, tStep);
570 
571  printf("Mean stress computed in Gnuplot export module: "); stress.printYourself();
572 
573  double time = 0.0;
574 
576  if ( ts != NULL ) {
577  time = ts->giveTargetTime();
578  }
579 
580  int bcIndex = iBC.giveNumber();
581 
582  std :: stringstream strMeanStress;
583  strMeanStress << "PrescribedGradientGnuplotMeanStress" << bcIndex << "Time" << time << ".dat";
584  std :: string nameMeanStress = strMeanStress.str();
585  std::vector<double> componentArray, stressArray;
586 
587  for(int i = 1; i <= stress.giveSize(); i++) {
588  componentArray.push_back(i);
589  stressArray.push_back(stress.at(i));
590  }
591 
592  XFEMDebugTools::WriteArrayToGnuplot(nameMeanStress, componentArray, stressArray);
593 
594  // Homogenized strain
595 
596  FloatArray grad;
597  iBC.giveGradientVoigt(grad);
598  outputGradient(iBC.giveNumber(), *iBC.giveDomain(), grad, tStep);
599 }
600 
602 {
603  FloatArray stress;
604  iBC.computeField(stress, tStep);
605 
606  printf("Mean stress computed in Gnuplot export module: "); stress.printYourself();
607  printf("sigXX: %.12e\n", stress(0) );
608 
609  double time = 0.0;
610 
612  if ( ts != NULL ) {
613  time = ts->giveTargetTime();
614  }
615 
616  int bcIndex = iBC.giveNumber();
617 
618  std :: stringstream strMeanStress;
619  strMeanStress << "PrescribedGradientGnuplotMeanStress" << bcIndex << "Time" << time << ".dat";
620  std :: string nameMeanStress = strMeanStress.str();
621  std::vector<double> componentArray, stressArray;
622 
623  for(int i = 1; i <= stress.giveSize(); i++) {
624  componentArray.push_back(i);
625  stressArray.push_back(stress.at(i));
626  }
627 
628  XFEMDebugTools::WriteArrayToGnuplot(nameMeanStress, componentArray, stressArray);
629 
630 
631  // Homogenized strain
632  FloatArray grad;
633  iBC.giveGradientVoigt(grad);
634  outputGradient(iBC.giveNumber(), *iBC.giveDomain(), grad, tStep);
635 
636 #if 0
637  FloatArray grad;
638  iBC.giveGradientVoigt(grad);
639  double timeFactor = iBC.giveTimeFunction()->evaluate(ts, VM_Total);
640  printf("timeFactor: %e\n", timeFactor );
641  grad.times(timeFactor);
642  printf("Mean grad computed in Gnuplot export module: "); grad.printYourself();
643 
644  std :: stringstream strMeanGrad;
645  strMeanGrad << "PrescribedGradientGnuplotMeanGrad" << bcIndex << "Time" << time << ".dat";
646  std :: string nameMeanGrad = strMeanGrad.str();
647  std::vector<double> componentArrayGrad, gradArray;
648 
649  for(int i = 1; i <= grad.giveSize(); i++) {
650  componentArrayGrad.push_back(i);
651  gradArray.push_back(grad.at(i));
652  }
653 
654  XFEMDebugTools::WriteArrayToGnuplot(nameMeanGrad, componentArrayGrad, gradArray);
655 #endif
656 
658 
659  // Traction node coordinates
660  std::vector< std::vector<FloatArray> > nodePointArray;
661  size_t numTracEl = iBC.giveNumberOfTractionElements();
662  for(size_t i = 0; i < numTracEl; i++) {
663 
664  std::vector<FloatArray> points;
665  FloatArray xS, xE;
666  iBC.giveTractionElCoord(i, xS, xE);
667  points.push_back(xS);
668  points.push_back(xE);
669 
670  nodePointArray.push_back(points);
671  }
672 
673  std :: stringstream strTractionNodes;
674  strTractionNodes << "TractionNodesGnuplotTime" << time << ".dat";
675  std :: string nameTractionNodes = strTractionNodes.str();
676 
677  WritePointsToGnuplot(nameTractionNodes, nodePointArray);
678 
679 
680 
681  // Traction element normal direction
682  std::vector< std::vector<FloatArray> > nodeNormalArray;
683  for(size_t i = 0; i < numTracEl; i++) {
684 
685  std::vector<FloatArray> points;
686  FloatArray n,t;
687  iBC.giveTractionElNormal(i, n,t);
688  points.push_back(n);
689  points.push_back(n);
690 
691  nodeNormalArray.push_back(points);
692  }
693 
694  std :: stringstream strTractionNodeNormals;
695  strTractionNodeNormals << "TractionNodeNormalsGnuplotTime" << time << ".dat";
696  std :: string nameTractionNodeNormals = strTractionNodeNormals.str();
697 
698  WritePointsToGnuplot(nameTractionNodeNormals, nodeNormalArray);
699 
700 
701 
702  // Traction (x,y)
703  std::vector< std::vector<FloatArray> > nodeTractionArray;
704  for(size_t i = 0; i < numTracEl; i++) {
705 
706  std::vector<FloatArray> tractions;
707  FloatArray tS, tE;
708 
709  iBC.giveTraction(i, tS, tE, VM_Total, tStep);
710 
711  tractions.push_back(tS);
712  tractions.push_back(tE);
713  nodeTractionArray.push_back(tractions);
714  }
715 
716  std :: stringstream strTractions;
717  strTractions << "TractionsGnuplotTime" << time << ".dat";
718  std :: string nameTractions = strTractions.str();
719 
720  WritePointsToGnuplot(nameTractions, nodeTractionArray);
721 
722 
723 
724  // Arc position along the boundary
725  std::vector< std::vector<FloatArray> > arcPosArray;
726  for(size_t i = 0; i < numTracEl; i++) {
727  std::vector<FloatArray> arcPos;
728  double xiS = 0.0, xiE = 0.0;
729  iBC.giveTractionElArcPos(i, xiS, xiE);
730  arcPos.push_back( FloatArray{xiS} );
731  arcPos.push_back( FloatArray{xiE} );
732 
733  arcPosArray.push_back(arcPos);
734  }
735 
736  std :: stringstream strArcPos;
737  strArcPos << "ArcPosGnuplotTime" << time << ".dat";
738  std :: string nameArcPos = strArcPos.str();
739 
740  WritePointsToGnuplot(nameArcPos, arcPosArray);
741 
742 
743  // Traction (normal, tangent)
744  std::vector< std::vector<FloatArray> > nodeTractionNTArray;
745  for(size_t i = 0; i < numTracEl; i++) {
746 
747  std::vector<FloatArray> tractions;
748  FloatArray tS, tE;
749 
750  iBC.giveTraction(i, tS, tE, VM_Total, tStep);
751  FloatArray n,t;
752  iBC.giveTractionElNormal(i, n, t);
753 
754 
755  double tSn = tS.dotProduct(n,2);
756  double tSt = tS.dotProduct(t,2);
757  tractions.push_back( {tSn ,tSt} );
758 
759  double tEn = tE.dotProduct(n,2);
760  double tEt = tE.dotProduct(t,2);
761  tractions.push_back( {tEn, tEt} );
762  nodeTractionNTArray.push_back(tractions);
763  }
764 
765  std :: stringstream strTractionsNT;
766  strTractionsNT << "TractionsNormalTangentGnuplotTime" << time << ".dat";
767  std :: string nameTractionsNT = strTractionsNT.str();
768 
769  WritePointsToGnuplot(nameTractionsNT, nodeTractionNTArray);
770 
771 
772 
773  // Boundary points and displacements
774  IntArray boundaries, bNodes;
775  iBC.giveBoundaries(boundaries);
776 
777  std::vector< std::vector<FloatArray> > bndNodes;
778 
779  for ( int pos = 1; pos <= boundaries.giveSize() / 2; ++pos ) {
780 
781  Element *e = iBC.giveDomain()->giveElement( boundaries.at(pos * 2 - 1) );
782  int boundary = boundaries.at(pos * 2);
783 
784  e->giveInterpolation()->boundaryGiveNodes(bNodes, boundary);
785 
786  std::vector<FloatArray> bndSegNodes;
787 
788  // Add the start and end nodes of the segment
789  DofManager *startNode = e->giveDofManager( bNodes[0] );
790  FloatArray xS = *(startNode->giveCoordinates());
791 
792  Dof *dSu = startNode->giveDofWithID(D_u);
793  double dU = dSu->giveUnknown(VM_Total, tStep);
794  xS.push_back(dU);
795 
796  Dof *dSv = startNode->giveDofWithID(D_v);
797  double dV = dSv->giveUnknown(VM_Total, tStep);
798  xS.push_back(dV);
799 
800  bndSegNodes.push_back(xS);
801 
802  DofManager *endNode = e->giveDofManager( bNodes[1] );
803  FloatArray xE = *(endNode->giveCoordinates());
804 
805  Dof *dEu = endNode->giveDofWithID(D_u);
806  dU = dEu->giveUnknown(VM_Total, tStep);
807  xE.push_back(dU);
808 
809  Dof *dEv = endNode->giveDofWithID(D_v);
810  dV = dEv->giveUnknown(VM_Total, tStep);
811  xE.push_back(dV);
812 
813  bndSegNodes.push_back(xE);
814 
815  bndNodes.push_back(bndSegNodes);
816  }
817 
818  std :: stringstream strBndNodes;
819  strBndNodes << "BndNodesGnuplotTime" << time << ".dat";
820  std :: string nameBndNodes = strBndNodes.str();
821 
822  WritePointsToGnuplot(nameBndNodes, bndNodes);
823 
824  }
825 }
826 
828 {
829  // Homogenized strain
830  double timeFactor = d.giveBc(bc)->giveTimeFunction()->evaluateAtTime(tStep->giveTargetTime());
831  printf("timeFactor: %e\n", timeFactor );
832  grad.times(timeFactor);
833  printf("Mean grad computed in Gnuplot export module: "); grad.printYourself();
834 
835  double time = tStep->giveTargetTime();
836 
837 
838  std :: stringstream strMeanGrad;
839  strMeanGrad << "PrescribedGradientGnuplotMeanGrad" << bc << "Time" << time << ".dat";
840  std :: string nameMeanGrad = strMeanGrad.str();
841  std::vector<double> componentArrayGrad, gradArray;
842 
843  for(int i = 1; i <= grad.giveSize(); i++) {
844  componentArrayGrad.push_back(i);
845  gradArray.push_back(grad.at(i));
846  }
847 
848  XFEMDebugTools::WriteArrayToGnuplot(nameMeanGrad, componentArrayGrad, gradArray);
849 }
850 
852 {
853  std::vector< std::vector<FloatArray> > pointArray;
854 
855  if(iDomain.giveNumberOfSpatialDimensions() == 2) {
856  // Write all element edges to gnuplot
857  for ( auto &el : iDomain.giveElements() ) {
858  int numEdges = el->giveNumberOfNodes();
859 
860 
861  for ( int edgeIndex = 1; edgeIndex <= numEdges; edgeIndex++ ) {
862  std::vector<FloatArray> points;
863 
864  IntArray bNodes;
865  el->giveInterpolation()->boundaryGiveNodes(bNodes, edgeIndex);
866 
867  int niLoc = bNodes.at(1);
868  const FloatArray &xS = *(el->giveNode(niLoc)->giveCoordinates() );
869  points.push_back(xS);
870 
871  int njLoc = bNodes.at( bNodes.giveSize() );
872  const FloatArray &xE = *(el->giveNode(njLoc)->giveCoordinates() );
873  points.push_back(xE);
874 
875  pointArray.push_back(points);
876  }
877 
878  }
879 
880 
881  double time = 0.0;
882 
884  if ( ts != NULL ) {
885  time = ts->giveTargetTime();
886  }
887 
888  std :: stringstream strMesh;
889  strMesh << "MeshGnuplotTime" << time << ".dat";
890  std :: string nameMesh = strMesh.str();
891 
892  WritePointsToGnuplot(nameMesh, pointArray);
893  }
894 }
895 
897 {
898  // Append node solution to history
899  FloatArray nodeSol;
900  iDMan.giveCompleteUnknownVector(nodeSol, VM_Total, tStep);
901 
902  mMonitorNodeDispHist.push_back(nodeSol);
903 
904 
905  // Write to file
906  std::vector< std::vector<FloatArray> > nodeDispHist;
907  nodeDispHist.push_back(mMonitorNodeDispHist);
908 
909  std :: string name = "MonitorNodeSolGnuplot.dat";
910  WritePointsToGnuplot(name, nodeDispHist);
911 }
912 
914 // printf("Exporting interface el.\n");
915 
916  std::vector<FloatArray> points, tractions, tractionsProj;
917 
918  int numEl = d.giveNumberOfElements();
919  for(int i = 1; i <= numEl; i++) {
920  Element *el = d.giveElement(i);
921 
922  StructuralInterfaceElement *intEl = dynamic_cast<StructuralInterfaceElement*>(el);
923 
924  if(intEl) {
925  printf("Found StructuralInterfaceElement.\n");
926 
928 
929  int numGP = ir->giveNumberOfIntegrationPoints();
930 // printf("numGP: %d\n", numGP );
931 
932  for(int gpInd = 0; gpInd < numGP; gpInd++ ) {
933  GaussPoint *gp = ir->getIntegrationPoint(gpInd);
934 // gp->giveGlobalCoordinates().printYourself();
935 
936  points.push_back( (gp->giveGlobalCoordinates()) );
937 
938  MaterialStatus *ms = static_cast< MaterialStatus * >( gp->giveMaterialStatus() );
940 
941  FloatArray traction = ims->giveTraction();
942 // printf("traction: "); traction.printYourself();
943  tractions.push_back(traction);
944 
945 
946  FloatArray tractionProj = ims->giveProjectedTraction();
947  tractionsProj.push_back(tractionProj);
948 
949  }
950 
951  }
952  }
953 
954  // Export x vs normal traction
955 
956  double time = 0.0;
957 
959  if ( ts != NULL ) {
960  time = ts->giveTargetTime();
961  }
962 
963  FILE * pFileX;
964  std :: stringstream strFileNameX;
965  strFileNameX << "NormalTractionVsXTime" << time << ".dat";
966  std :: string nameStringX = strFileNameX.str();
967 
968  pFileX = fopen ( nameStringX.c_str() , "wb" );
969 
970  fprintf(pFileX, "#x tn\n");
971  for ( size_t j = 0; j < points.size(); j++ ) {
972  fprintf(pFileX, "%e %e\n", points[j][0], tractions[j].at(3) );
973  }
974 
975  fclose(pFileX);
976 
977 
978 // FILE * pFileXProj;
979 // std :: stringstream strFileNameXProj;
980 // strFileNameXProj << "NormalTractionProjVsXTime" << time << ".dat";
981 // std :: string nameStringXProj = strFileNameXProj.str();
982 //
983 // pFileXProj = fopen ( nameStringXProj.c_str() , "wb" );
984 //
985 // fprintf(pFileXProj, "#x tn\n");
986 // for ( size_t j = 0; j < points.size(); j++ ) {
988 // fprintf(pFileXProj, "%e %e\n", points[j][0], tractionsProj[j].at(3) );
989 // }
990 //
991 // fclose(pFileXProj);
992 
993 
994  // Export x vs shear traction
995 
996  FILE * pFileXshear;
997  std :: stringstream strFileNameXshear;
998  strFileNameXshear << "ShearTractionVsXTime" << time << ".dat";
999  std :: string nameStringXshear = strFileNameXshear.str();
1000 
1001  pFileXshear = fopen ( nameStringXshear.c_str() , "wb" );
1002 
1003  fprintf(pFileXshear, "#x tn\n");
1004  for ( size_t j = 0; j < points.size(); j++ ) {
1005  fprintf(pFileXshear, "%e %e\n", points[j][0], tractions[j].at(1) );
1006  }
1007 
1008  fclose(pFileXshear);
1009 
1010 
1011 
1012 
1013  // Export y vs normal traction
1014  FILE * pFileY;
1015  std :: stringstream strFileNameY;
1016  strFileNameY << "NormalTractionVsYTime" << time << ".dat";
1017  std :: string nameStringY = strFileNameY.str();
1018 
1019  pFileY = fopen ( nameStringY.c_str() , "wb" );
1020 
1021  fprintf(pFileY, "#y tn\n");
1022  for ( size_t j = 0; j < points.size(); j++ ) {
1023  fprintf(pFileY, "%e %e\n", points[j][1], tractions[j].at(3) );
1024  }
1025 
1026  fclose(pFileY);
1027 
1028 
1029 
1030 
1031  // Export y vs shear traction
1032  FILE * pFileYshear;
1033  std :: stringstream strFileNameYshear;
1034  strFileNameYshear << "ShearTractionVsYTime" << time << ".dat";
1035  std :: string nameStringYshear = strFileNameYshear.str();
1036 
1037  pFileYshear = fopen ( nameStringYshear.c_str() , "wb" );
1038 
1039  fprintf(pFileYshear, "#y tn\n");
1040  for ( size_t j = 0; j < points.size(); j++ ) {
1041  fprintf(pFileYshear, "%e %e\n", points[j][1], tractions[j].at(1) );
1042  }
1043 
1044  fclose(pFileYshear);
1045 
1046 }
1047 
1048 
1049 void GnuplotExportModule :: WritePointsToGnuplot(const std :: string &iName, const std :: vector< std::vector<FloatArray> > &iPoints)
1050 {
1051  std :: ofstream file;
1052  file.open( iName.data() );
1053 
1054  // Set some output options
1055  file << std :: scientific;
1056 
1057  file << "# x y\n";
1058 
1059  for(auto posVec: iPoints) {
1060  for(auto pos: posVec) {
1061 
1062  for(int i = 0; i < pos.giveSize(); i++) {
1063  file << pos[i] << " ";
1064  }
1065  file << "\n";
1066 
1067 // file << pos[0] << " " << pos[1] << "\n";
1068  }
1069  file << "\n";
1070  }
1071 
1072  file.close();
1073 
1074 }
1075 
1076 
1077 } // end namespace oofem
virtual void computeField(FloatArray &sigma, TimeStep *tStep)
Computes the homogenized, macroscopic, field (stress).
bool testTimeStepOutput(TimeStep *tStep)
Tests if given time step output is required.
Definition: exportmodule.C:148
Prescribes or where are primary unknowns for the subscale.
std::vector< double > mTimeHist
Abstract class representing entity, which is included in the FE model using one (or more) global func...
void outputBoundaryCondition(PrescribedGradient &iBC, TimeStep *tStep)
Boundary condition output.
EnrichmentFront * giveEnrichmentFrontStart()
Class and object Domain.
Definition: domain.h:115
void giveSubPolygon(std::vector< FloatArray > &oPoints, const double &iXiStart, const double &iXiEnd) const
virtual IntegrationRule * giveDefaultIntegrationRulePtr()
Access method for default integration rule.
Definition: element.h:822
const std::vector< double > & giveCohesiveZoneArcPositions() const
Definition: crack.h:69
void giveTractionElCoord(size_t iElInd, FloatArray &oStartCoord, FloatArray &oEndCoord) const
virtual double give(Dof *dof, ValueModeType mode, TimeStep *tStep)
Returns the value of a prescribed unknown, respecting requested mode for given time.
void outputXFEM(EnrichmentItem &iEI, TimeStep *tStep)
XFEM output.
int giveNumberOfBoundaryConditions() const
Returns number of boundary conditions in domain.
Definition: domain.h:440
Imposes a prescribed gradient weakly on the boundary with a Neumann boundary condition.
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
void outputXFEMGeometry(const std::vector< std::vector< FloatArray > > &iEnrItemPoints)
TipInfo gathers useful information about a crack tip, like its position and tangent direction...
Definition: tipinfo.h:24
virtual void computeTangentialSignDist(double &oDist, const FloatArray &iPoint, double &oMinDistArcPos) const =0
virtual FloatArray * giveCoordinates()
Definition: dofmanager.h:382
virtual double giveUnknown(ValueModeType mode, TimeStep *tStep)=0
The key method of class Dof.
virtual bool hasField(InputFieldType id)=0
Returns true if record contains field identified by idString keyword.
virtual void initialize()
Definition: exportmodule.C:86
void outputGradient(int bc, Domain &d, FloatArray &grad, TimeStep *tStep)
Help functions.
double giveTargetTime()
Returns target time.
Definition: timestep.h:146
const FloatArray & giveFirstPKTraction() const
Returns the const pointer to receiver&#39;s first Piola-Kirchhoff traction vector.
Abstract base class for all finite elements.
Definition: element.h:145
void giveTraction(size_t iElInd, FloatArray &oStartTraction, FloatArray &oEndTraction, ValueModeType mode, TimeStep *tStep)
Base class for dof managers.
Definition: dofmanager.h:113
int giveNumberOfElements() const
Returns number of elements in domain.
Definition: domain.h:434
double evaluate(TimeStep *tStep, ValueModeType mode)
Returns the value of load time function at given time.
Definition: function.C:55
int giveNumberOfSpatialDimensions()
Returns number of spatial dimensions.
Definition: domain.C:1067
Represents export output module - a base class for all output modules.
Definition: exportmodule.h:71
Class implementing an array of integers.
Definition: intarray.h:61
int & at(int i)
Coefficient access function.
Definition: intarray.h:103
Abstract representation of Geometry.
Definition: geometry.h:81
XfemManager * giveXfemManager()
Definition: domain.C:375
virtual void terminate()
Terminates the receiver.
virtual FEInterpolation * giveInterpolation() const
Definition: element.h:629
#define _IFT_GnuplotExportModule_xfem
Abstract base class representing integration rule.
void giveTractionElNormal(size_t iElInd, FloatArray &oNormal, FloatArray &oTangent) const
GeneralBoundaryCondition * giveBc(int n)
Service for accessing particular domain bc.
Definition: domain.C:243
#define _IFT_GnuplotExportModule_interface_el
virtual void computeField(FloatArray &sigma, TimeStep *tStep)
Computes the homogenized, macroscopic, field (stress).
void outputInterfaceEl(Domain &d, TimeStep *tStep)
Custom output for interface elements.
int giveNumberOfEnrichmentItems() const
Definition: xfemmanager.h:185
#define _IFT_GnuplotExportModule_BoundaryConditionsExtra
void buildReactionTable(IntArray &restrDofMans, IntArray &restrDofs, IntArray &eqn, TimeStep *tStep, int di)
Builds the reaction force table.
#define max
Crack.
Definition: crack.h:54
virtual int giveBcId()=0
Returns the id of associated boundary condition, if there is any.
Class implementing Dirichlet boundary condition on DOF (primary boundary condition).
void giveBoundaries(IntArray &oBoundaries)
Element * giveElement(int n)
Service for accessing particular domain fe element.
Definition: domain.C:160
double dotProduct(const FloatArray &x) const
Computes the dot product (or inner product) of receiver and argument.
Definition: floatarray.C:463
#define OOFEM_ERROR(...)
Definition: error.h:61
#define _IFT_GnuplotExportModule_BoundaryConditions
EngngModel * emodel
Problem pointer.
Definition: exportmodule.h:77
static void WriteArrayToGnuplot(const std::string &iName, const std::vector< double > &iX, const std::vector< double > &iY)
const FloatArray & giveNormal() const
Class EnrichmentFront: describes the edge or tip of an XFEM enrichment.
static void WritePointsToGnuplot(const std::string &iName, const std::vector< std::vector< FloatArray > > &iPoints)
This class implements a structural interface material status information.
#define _IFT_GnuplotExportModule_monitornode
void giveCompleteUnknownVector(FloatArray &answer, ValueModeType mode, TimeStep *tStep)
Assembles the complete unknown vector in node.
Definition: dofmanager.C:737
Imposes a prescribed gradient weakly on the boundary with an independent traction discretization...
const std::vector< GaussPoint * > & giveCohesiveZoneGaussPoints() const
Definition: crack.h:68
const FloatArray & giveJump() const
Returns the const pointer to receiver&#39;s jump.
const FloatArray & giveGlobalCoordinates()
Definition: gausspoint.h:160
Abstract base class representing a material status information.
Definition: matstatus.h:84
bool hasXfemManager()
Definition: domain.C:386
virtual void boundaryGiveNodes(IntArray &answer, int boundary)=0
Gives the boundary nodes for requested boundary number.
Class representing vector of real numbers.
Definition: floatarray.h:82
const FloatArray & giveTraction() const
Returns the const pointer to receiver&#39;s traction vector.
virtual void computeField(FloatArray &sigma, TimeStep *tStep)
Computes the homogenized, macroscopic, field (stress).
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
This class manages the xfem part.
Definition: xfemmanager.h:109
IRResultType
Type defining the return values of InputRecord reading operations.
Definition: irresulttype.h:47
GaussPoint * getIntegrationPoint(int n)
Access particular integration point of receiver.
const FloatArray & giveStressVector() const
Returns the const pointer to receiver&#39;s stress vector.
Definition: structuralms.h:107
std::unique_ptr< MaterialForceEvaluator > mpMatForceEvaluator
Evaluator for material forces.
#define _IFT_GnuplotExportModule_cracklength
IntegrationPointStatus * giveMaterialStatus()
Returns reference to associated material status (NULL if not defined).
Definition: gausspoint.h:205
std::vector< FloatArray > mMonitorNodeDispHist
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: exportmodule.C:58
#define _IFT_GnuplotExportModule_ReactionForces
const FloatArray & giveProjectedTraction() const
Returns the projected traction.
Class representing the general Input Record.
Definition: inputrecord.h:101
virtual void printYourself() const
Print receiver on stdout.
Definition: floatarray.C:747
virtual void doOutput(TimeStep *tStep, bool forcedOutput=false)
Writes the output.
BasicGeometry * giveGeometry()
Dof * giveDofWithID(int dofID) const
Returns DOF with given dofID; issues error if not present.
Definition: dofmanager.C:119
int giveNumberOfIntegrationPoints() const
Returns number of integration points of receiver.
double computeLength()
Definition: crack.C:149
void giveTractionElArcPos(size_t iElInd, double &oXiStart, double &oXiEnd) const
void times(double s)
Multiplies receiver with scalar.
Definition: floatarray.C:818
EnrichmentItem * giveEnrichmentItem(int n)
Definition: xfemmanager.h:184
std::vector< std::vector< FloatArray > > mReactionForceHistory
Stores the sum of reaction forces for each BC.
#define new
void outputReactionForces(TimeStep *tStep)
void giveGradientVoigt(FloatArray &oGradient) const
Gives back the applied gradient in Voigt form.
Domain * giveDomain() const
Definition: femcmpnn.h:100
This class implements extension of EngngModel for structural models.
#define _IFT_GnuplotExportModule_materialforceradii
const TipInfo & giveTipInfo() const
#define _IFT_GnuplotExportModule_mesh
void push_back(const double &iVal)
Add one element.
Definition: floatarray.h:118
std::vector< std::unique_ptr< Element > > & giveElements()
Definition: domain.h:279
void computeReaction(FloatArray &answer, TimeStep *tStep, int di)
Computes reaction forces.
Abstract base class representing the "problem" under consideration.
Definition: engngm.h:181
Abstract base class for all structural interface elements.
virtual TimeStep * giveCurrentStep(bool force=false)
Returns current time step.
Definition: engngm.h:683
int giveSize() const
Definition: intarray.h:203
Evaluates material forces.
int giveSize() const
Returns the size of receiver.
Definition: floatarray.h:218
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
std::vector< std::vector< double > > mDispHist
the oofem namespace is to define a context or scope in which all oofem names are defined.
DofManager * giveDofManager(int i) const
Definition: element.C:514
Domain * giveDomain(int n)
Service for accessing particular problem domain.
Definition: engngm.C:1720
Abstract class Dof represents Degree Of Freedom in finite element mesh.
Definition: dof.h:93
virtual void callGnuplotExportModule(GnuplotExportModule &iExpMod, TimeStep *tStep)
int giveNumber() const
Definition: femcmpnn.h:107
DofManager * giveDofManager(int n)
Service for accessing particular domain dof manager.
Definition: domain.C:314
void outputMesh(Domain &iDomain)
Mesh output.
void outputNodeDisp(DofManager &iDMan, TimeStep *tStep)
Monitor node output.
virtual double evaluateAtTime(double t)
Returns the value of the function at given time.
Definition: function.C:76
(Under development) The Gnuplot export module enables OOFEM to export some data in a format that can ...
Class representing integration point in finite element program.
Definition: gausspoint.h:93
IRResultType giveOptionalField(int &answer, InputFieldType id)
Reads the integer field value.
Definition: inputrecord.C:64
std::unordered_map< int, std::vector< double > > mCrackLengthHist
Store time history of crack lengths.
Class representing solution step.
Definition: timestep.h:80
REGISTER_ExportModule(ErrorCheckingExportModule)
EnrichmentFront * giveEnrichmentFrontEnd()
EnrichmentItem with geometry described by BasicGeometry.

This page is part of the OOFEM documentation. Copyright (c) 2011 Borek Patzak
Project e-mail: info@oofem.org
Generated at Tue Jan 2 2018 20:07:29 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011