4 * Created on: 7 fevrier. 2012
11 #include "MEDFileMesh.hxx"
12 #include "MEDLoader.hxx"
13 #include "MEDCouplingUMesh.hxx"
14 #include "MEDCouplingFieldDouble.hxx"
16 #include "CdmathException.hxx"
20 using namespace MEDCoupling;
24 //----------------------------------------------------------------------
25 Field::Field( EntityType typeField )
26 //----------------------------------------------------------------------
30 _numberOfComponents=0;
33 //----------------------------------------------------------------------
35 //----------------------------------------------------------------------
37 //std::cerr << "dtor Field, _field = " <<_field << std::endl;
38 //if (_field) _field->decrRef();
42 Field::Field(const std::string fieldName, EntityType type, const Mesh& mesh, int numberOfComponents, double time)
47 _numberOfComponents=numberOfComponents;
51 buildFieldMemoryStructure();
53 void Field::buildFieldMemoryStructure()
55 MEDCouplingUMesh* mu=_mesh.getMEDCouplingMesh()->buildUnstructured();
56 DataArrayDouble *array=DataArrayDouble::New();
57 if (_typeField==CELLS)
59 _field=MEDCouplingFieldDouble::New(ON_CELLS);
60 array->alloc(_mesh.getNumberOfCells(),_numberOfComponents);
62 }else if(_typeField==NODES)
64 _field=MEDCouplingFieldDouble::New(ON_NODES);
65 array->alloc(_mesh.getNumberOfNodes(),_numberOfComponents);
67 }else if(_typeField==FACES)
69 _field=MEDCouplingFieldDouble::New(ON_CELLS);
70 array->alloc(_mesh.getNumberOfFaces(),_numberOfComponents);
71 DataArrayInt *desc=DataArrayInt::New();
72 DataArrayInt *descI=DataArrayInt::New();
73 DataArrayInt *revDesc=DataArrayInt::New();
74 DataArrayInt *revDescI=DataArrayInt::New();
75 MEDCouplingUMesh *m3=mu->buildDescendingConnectivity(desc,descI,revDesc,revDescI);
83 throw CdmathException("Type of Field::Field() is not compatible");
85 _field->setName(_fieldName.c_str()) ;
86 _field->setArray(array);
87 _field->setTime(_time,0,0);
92 Field::Field( const std::string filename, EntityType type,
93 const std::string & fieldName,
94 int iteration, int order, int meshLevel,
95 int numberOfComponents, double time)
98 _mesh=Mesh(filename + ".med", meshLevel);
100 _numberOfComponents=numberOfComponents;
102 _fieldName=fieldName;
104 readFieldMed(filename, type, fieldName, iteration, order);
107 Field::Field(const std::string meshFileName, EntityType type, const std::vector<double> Vconstant,
108 const std::string & fieldName, int meshLevel, double time )
111 _mesh=Mesh(meshFileName + ".med", meshLevel);
113 _numberOfComponents=Vconstant.size();
115 _fieldName=fieldName;
117 buildFieldMemoryStructure();
119 int nbelem=_field->getNumberOfTuples();
120 int nbcomp=_field->getNumberOfComponents() ;
122 for (int ielem=0 ; ielem<nbelem; ielem++)
123 for (int jcomp=0 ; jcomp<nbcomp ; jcomp++)
124 _field->getArray()->getPointer()[jcomp+ielem*_field->getNumberOfComponents()]=Vconstant[jcomp];
126 Field::Field(const Mesh& M, EntityType type, const Vector Vconstant, const std::string & fieldName, double time)
131 _numberOfComponents=Vconstant.size();
133 _fieldName=fieldName;
135 buildFieldMemoryStructure();
137 int nbelem=_field->getNumberOfTuples();
138 int nbcomp=_field->getNumberOfComponents() ;
140 for (int ielem=0 ; ielem<nbelem; ielem++)
141 for (int jcomp=0 ; jcomp<nbcomp ; jcomp++)
142 _field->getArray()->getPointer()[jcomp+ielem*_field->getNumberOfComponents()]=Vconstant[jcomp];
144 Field::Field(const Mesh& M, EntityType type, const vector<double> Vconstant, const std::string & fieldName, double time)
149 _numberOfComponents=Vconstant.size();
151 _fieldName=fieldName;
153 buildFieldMemoryStructure();
155 int nbelem=_field->getNumberOfTuples();
156 int nbcomp=_field->getNumberOfComponents() ;
158 for (int ielem=0 ; ielem<nbelem; ielem++)
159 for (int jcomp=0 ; jcomp<nbcomp ; jcomp++)
160 _field->getArray()->getPointer()[jcomp+ielem*_field->getNumberOfComponents()]=Vconstant[jcomp];
162 Field::Field( int nDim, const vector<double> Vconstant, EntityType type,
163 double xmin, double xmax, int nx, string leftSide, string rightSide,
164 double ymin, double ymax, int ny, string backSide, string frontSide,
165 double zmin, double zmax, int nz, string bottomSide, string topSide,
166 const std::string & fieldName, double time,double epsilon)
170 _numberOfComponents=Vconstant.size();
172 _fieldName=fieldName;
176 _mesh=Mesh(xmin,xmax,nx);
179 _mesh=Mesh(xmin,xmax,nx,ymin,ymax,ny);
181 _mesh=Mesh(xmin,xmax,nx,ymin,ymax,ny,zmin,zmax,nz);
183 cout<<"Field::Field: Space dimension nDim should be between 1 and 3"<<endl;
184 throw CdmathException("Space dimension nDim should be between 1 and 3");
187 _mesh.setGroupAtPlan(xmax,0,epsilon,rightSide);
188 _mesh.setGroupAtPlan(xmin,0,epsilon,leftSide);
190 _mesh.setGroupAtPlan(ymax,1,epsilon,frontSide);
191 _mesh.setGroupAtPlan(ymin,1,epsilon,backSide);
194 _mesh.setGroupAtPlan(zmax,2,epsilon,topSide);
195 _mesh.setGroupAtPlan(zmin,2,epsilon,bottomSide);
199 buildFieldMemoryStructure();
201 int nbelem=_field->getNumberOfTuples();
202 int nbcomp=_field->getNumberOfComponents() ;
204 for (int ielem=0 ; ielem<nbelem; ielem++)
205 for (int jcomp=0 ; jcomp<nbcomp ; jcomp++)
206 _field->getArray()->getPointer()[jcomp+ielem*_field->getNumberOfComponents()]=Vconstant[jcomp];
208 Field::Field(const Mesh M, const Vector VV_Left, const Vector VV_Right, double disc_pos,
209 EntityType type, int direction, const std::string & fieldName, double time)
211 if (VV_Right.getNumberOfRows()!=VV_Left.getNumberOfRows())
212 throw CdmathException( "Field::Field: Vectors VV_Left and VV_Right have different sizes");
217 _numberOfComponents=VV_Left.getNumberOfRows();
219 _fieldName=fieldName;
222 buildFieldMemoryStructure();
224 int nbelem=_field->getNumberOfTuples();
225 int nbcomp=_field->getNumberOfComponents() ;
226 double component_value;
228 for (int j = 0; j < nbelem; j++) {
230 component_value=M.getCell(j).x();
231 else if(direction==1)
232 component_value=M.getCell(j).y();
233 else if(direction==2)
234 component_value=M.getCell(j).z();
236 throw CdmathException( "Field::Field: direction should be an integer between 0 and 2");
238 for (int i=0; i< nbcomp; i++)
239 if (component_value< disc_pos )
240 _field->getArray()->getPointer()[j+i*_field->getNumberOfComponents()] = VV_Left[i];
242 _field->getArray()->getPointer()[j+i*_field->getNumberOfComponents()] = VV_Right[i];
245 Field::Field( int nDim, const vector<double> VV_Left, vector<double> VV_Right,
246 double xstep, EntityType type,
247 double xmin, double xmax, int nx, string leftSide, string rightSide,
248 double ymin, double ymax, int ny, string backSide, string frontSide,
249 double zmin, double zmax, int nz, string bottomSide, string topSide,
250 int direction, const std::string & fieldName, double time, double epsilon)
252 if (VV_Right.size()!=VV_Left.size())
253 throw CdmathException( "Field::Field: Vectors VV_Left and VV_Right have different sizes");
256 M=Mesh(xmin,xmax,nx);
258 M=Mesh(xmin,xmax,nx,ymin,ymax,ny);
260 M=Mesh(xmin,xmax,nx,ymin,ymax,ny,zmin,zmax,nz);
262 throw CdmathException("Field::Field : Space dimension nDim should be between 1 and 3");
264 M.setGroupAtPlan(xmax,0,epsilon,rightSide);
265 M.setGroupAtPlan(xmin,0,epsilon,leftSide);
267 M.setGroupAtPlan(ymax,1,epsilon,frontSide);
268 M.setGroupAtPlan(ymin,1,epsilon,backSide);
271 M.setGroupAtPlan(zmax,2,epsilon,topSide);
272 M.setGroupAtPlan(zmin,2,epsilon,bottomSide);
274 Vector V_Left(VV_Left.size()), V_Right(VV_Right.size());
275 for(int i=0;i<VV_Left.size(); i++){
276 V_Left(i)=VV_Left[i];
277 V_Right(i)=VV_Right[i];
279 Field(M, V_Left, V_Right, xstep, type, direction, fieldName, time);
282 Field::Field(const Mesh M, const Vector Vin, const Vector Vout, double radius,
283 const Vector Center, EntityType type, const std::string & fieldName, double time)
285 if((Center.size()!=M.getSpaceDimension()) || (Vout.size() != Vin.size()) )
287 cout<< "Vout.size()= "<<Vout.size() << ", Vin.size()= "<<Vin.size()<<", Center.size()="<<Center.size()<<", M.getSpaceDim= "<< M.getSpaceDimension()<<endl;
288 throw CdmathException("Field::Field : Vector size error");
294 _numberOfComponents=Vout.size();
296 _fieldName=fieldName;
299 buildFieldMemoryStructure();
301 int nbelem=_field->getNumberOfTuples();
302 int nbcomp=_field->getNumberOfComponents() ;
304 int spaceDim=M.getSpaceDimension();
305 Vector currentPoint(spaceDim);
307 for(int i=0;i<nbelem;i++)
309 currentPoint(0)=M.getCell(i).x();
312 currentPoint(1)=M.getCell(i).y();
314 currentPoint(2)=M.getCell(i).z();
316 if((currentPoint-Center).norm()<radius)
317 for(int j=0;j<nbcomp;j++)
318 _field->getArray()->getPointer()[j+i*_field->getNumberOfComponents()]=Vin[j];
320 for(int j=0;j<nbcomp;j++)
321 _field->getArray()->getPointer()[j+i*_field->getNumberOfComponents()]=Vout[j];
325 MEDCoupling::DataArrayDouble * Field::getArray(){
326 return _field->getArray();
330 Field::readFieldMed( const std::string & fileNameRadical,
332 const std::string & fieldName,
337 * Reads the file fileNameRadical.med and creates a Field from it.
339 std::string completeFileName = fileNameRadical + ".med";
340 std::vector<std::string> fieldNames = MEDCoupling::GetAllFieldNames(completeFileName);
342 std::string attributedFieldName;
345 // Get the name of the right field that we will attribute to the Field.
346 if (fieldName == "") {
347 if (fieldNames.size() > 0)
348 attributedFieldName = fieldNames[0];
350 std::ostringstream message;
351 message << "No field in file " << completeFileName;
352 throw CdmathException(message.str().c_str());
356 for (; iField < fieldNames.size(); iField++)
357 if (fieldName == fieldNames[iField]) break;
359 if (iField < fieldNames.size())
360 attributedFieldName = fieldName;
362 std::ostringstream message;
363 message << "No field named " << fieldName << " in file " << completeFileName;
364 throw CdmathException(message.str().c_str());
368 // Get the name of the right mesh that we will attribute to the Field.
369 std::vector<std::string> meshNames
370 = MEDCoupling::GetMeshNamesOnField(completeFileName, attributedFieldName);
371 if (meshNames.size() == 0) {
372 std::ostringstream message;
373 message << "No mesh associated to " << fieldName
374 << " in file " << completeFileName;
375 throw CdmathException(message.str().c_str());
377 std::string attributedMeshName = meshNames[0];
380 MEDCoupling::TypeOfField medFieldType[3] = { ON_CELLS, ON_NODES, ON_CELLS };
383 _field = dynamic_cast< MEDCoupling::MEDCouplingFieldDouble * > (
384 MEDCoupling::ReadFieldCell( completeFileName,
385 attributedMeshName, 0,
386 attributedFieldName, iteration, order) );
389 _field = dynamic_cast< MEDCoupling::MEDCouplingFieldDouble * > (
390 MEDCoupling::ReadFieldNode( completeFileName,
391 attributedMeshName, 0,
392 attributedFieldName, iteration, order) );
395 _field = dynamic_cast< MEDCoupling::MEDCouplingFieldDouble * > (
396 MEDCoupling::ReadFieldCell( completeFileName,
397 attributedMeshName, -1,
398 attributedFieldName, iteration, order) );
405 Field::getNormEuclidean() const
407 DoubleTab norm(getNumberOfElements(),_field->magnitude()->getArray()->getConstPointer());
414 if( getNumberOfComponents() !=1)
415 throw CdmathException("double Field::max() : field should have a single component in order to extract maximum value");
418 for(int i=0; i<getNumberOfElements() ; i++)
419 if( result < (*this)(i,0))
420 result = (*this)(i,0);
428 if( getNumberOfComponents() !=1)
429 throw CdmathException("double Field::min() : field should have a single component in order to extract minimum value");
432 for(int i=0; i<getNumberOfElements() ; i++)
433 if( result > (*this)(i,0))
434 result = (*this)(i,0);
440 Field::getInfoOnComponent(int icomp) const
442 return _field->getArray()->getInfoOnComponent(icomp);
446 Field::setInfoOnComponent(int icomp, string nameCompo)
448 _field.retn()->getArray()->setInfoOnComponent(icomp,nameCompo);
451 //----------------------------------------------------------------------
452 Field::Field( const Field & f )
453 //----------------------------------------------------------------------
456 MCAuto<MEDCouplingFieldDouble> f1=f.getField()->deepCopy();
458 _typeField=f.getTypeOfField();
461 //----------------------------------------------------------------------
462 MCAuto<MEDCouplingFieldDouble>
463 Field::getField ( void ) const
464 //----------------------------------------------------------------------
469 //----------------------------------------------------------------------
471 Field::setFieldByMEDCouplingFieldDouble ( const MEDCouplingFieldDouble* field )
472 //----------------------------------------------------------------------
474 MCAuto<MEDCouplingFieldDouble> ff=field->deepCopy();
478 //----------------------------------------------------------------------
480 Field::setFieldByDataArrayDouble ( const DataArrayDouble* array )
481 //----------------------------------------------------------------------
483 _field->setArray(const_cast<DataArrayDouble*>(array));
486 //----------------------------------------------------------------------
488 Field::integral ( ) const
489 //----------------------------------------------------------------------
491 int nbComp=_field->getNumberOfComponents();
492 double res[nbComp];//Pointer containing the inegral of each component
493 _field->integral(false,res);
494 Vector result(nbComp);//Vector containing the inegral of each component
496 for(int i=0; i<nbComp ; i++)
502 //----------------------------------------------------------------------
504 Field::integral ( int compId ) const
505 //----------------------------------------------------------------------
507 return _field->integral(compId, false);
510 //----------------------------------------------------------------------
512 Field::normL1 ( ) const
513 //----------------------------------------------------------------------
515 int nbComp=_field->getNumberOfComponents();
516 double res[nbComp];//Pointer containing the L1 norm of each component
518 Vector result(nbComp);//Vector containing the L1 norm of each component
520 for(int i=0; i<nbComp ; i++)
526 //----------------------------------------------------------------------
528 Field::normL2 ( ) const
529 //----------------------------------------------------------------------
531 int nbComp=_field->getNumberOfComponents();
532 double res[nbComp];//Pointer containing the L2 norm of each component
534 Vector result(nbComp);//Vector containing the L2 norm of each component
536 for(int i=0; i<nbComp ; i++)
542 //----------------------------------------------------------------------
544 Field::normMax ( ) const
545 //----------------------------------------------------------------------
547 int nbComp=_field->getNumberOfComponents();
548 int nbElems=getNumberOfElements();
550 Vector result(nbComp);//Vector containing the Linfinity norm of each component
552 for(int i=0; i<nbElems ; i++)
553 for(int j=0; j<nbComp ; j++)
554 if(fabs((*this)(i,j))>result(j))
555 result(j)=fabs((*this)(i,j));
560 //----------------------------------------------------------------------
562 Field::operator() ( int ielem )
563 //----------------------------------------------------------------------
565 if(ielem>_field->getNumberOfTuples() || ielem<0)
566 throw CdmathException("double& Field::operator(ielem) : ielem>number of values !");
567 return _field->getArray()->getPointer()[ielem*_field->getNumberOfComponents()];
570 //----------------------------------------------------------------------
572 Field::operator[] ( int ielem )
573 //----------------------------------------------------------------------
575 if(ielem>_field->getNumberOfTuples() || ielem<0)
576 throw CdmathException("double& Field::operator[ielem] : ielem>number of values !");
577 return _field->getArray()->getPointer()[ielem*_field->getNumberOfComponents()];
580 //----------------------------------------------------------------------
582 Field::operator() ( int ielem ) const
583 //----------------------------------------------------------------------
585 if(ielem>_field->getNumberOfTuples() || ielem<0)
586 throw CdmathException("double Field::operator(ielem) : ielem>number of values !");
587 return _field->getArray()->getConstPointer()[ielem*_field->getNumberOfComponents()];
590 //----------------------------------------------------------------------
592 Field::operator[] ( int ielem ) const
593 //----------------------------------------------------------------------
595 if(ielem>_field->getNumberOfTuples() || ielem<0)
596 throw CdmathException("double Field::operator[ielem] : ielem>number of values !");
597 return _field->getArray()->getConstPointer()[ielem*_field->getNumberOfComponents()];
600 //----------------------------------------------------------------------
602 Field::operator() ( int ielem, int jcomp )
603 //----------------------------------------------------------------------
605 if(ielem>_field->getNumberOfTuples() || jcomp>_field->getNumberOfComponents() || ielem<0 || jcomp<0)
606 throw CdmathException("double& Field::operator( int ielem, int jcomp ) : ielem>number of values or jcomp>number of components !");
607 return _field->getArray()->getPointer()[jcomp+ielem*_field->getNumberOfComponents()];
610 //----------------------------------------------------------------------
612 Field::operator() ( int ielem, int jcomp ) const
613 //----------------------------------------------------------------------
615 if(ielem>_field->getNumberOfTuples() || jcomp>_field->getNumberOfComponents() || ielem<0 || jcomp<0)
616 throw CdmathException("double Field::operator( int ielem, int jcomp ) : ielem>number of values or jcomp>number of components !");
617 return _field->getArray()->getConstPointer()[jcomp+ielem*_field->getNumberOfComponents()];
620 //----------------------------------------------------------------------
622 Field::setTime ( double time, int iter )
623 //----------------------------------------------------------------------
625 _field->setTime(time,iter,0.0);
627 //----------------------------------------------------------------------
629 Field::getTime ( void ) const
630 //----------------------------------------------------------------------
633 return _field->getTime(a,b);
636 //----------------------------------------------------------------------
638 Field::getNumberOfElements ( void ) const
639 //----------------------------------------------------------------------
641 return _field->getNumberOfTuples() ;
645 Field::getSpaceDimension( void ) const
647 return _mesh.getSpaceDimension() ;
650 //----------------------------------------------------------------------
652 Field::getNumberOfComponents ( void ) const
653 //----------------------------------------------------------------------
655 return _field->getNumberOfComponents() ;
658 //----------------------------------------------------------------------
660 Field::getValues ( void ) const
661 //----------------------------------------------------------------------
663 return _field->getArray()->getConstPointer() ;
666 //----------------------------------------------------------------------
668 Field::getName ( void ) const
669 //----------------------------------------------------------------------
671 return _field->getName() ;
674 //----------------------------------------------------------------------
676 Field::getMesh ( void ) const
677 //----------------------------------------------------------------------
682 //----------------------------------------------------------------------
684 Field::getTypeOfField ( void ) const
685 //----------------------------------------------------------------------
690 //----------------------------------------------------------------------
692 Field::setName ( const string fieldName )
693 //----------------------------------------------------------------------
695 _field->setName(fieldName.c_str()) ;
699 //----------------------------------------------------------------------
701 Field::operator+ ( const Field& f ) const
702 //----------------------------------------------------------------------
704 //if(f.getMesh().getMEDCouplingMesh() != _mesh.getMEDCouplingMesh())
705 //throw CdmathException("Field::operator+ : Field addition requires identical meshes");
706 _mesh.getMEDCouplingMesh()->checkFastEquivalWith(f.getMesh().getMEDCouplingMesh(),1e-6);
707 if(f.getTypeOfField() != getTypeOfField())
708 throw CdmathException("Field::operator+ : Field addition requires identical field types (CELLS, NODES or FACES");
709 if(f.getNumberOfComponents() != getNumberOfComponents())
710 throw CdmathException("Field::operator+ : Field addition requires identical number of components");
712 Field fres(getName(),f.getTypeOfField(),f.getMesh(),f.getNumberOfComponents(),f.getTime());
713 int nbComp=f.getNumberOfComponents();
714 int nbElem=f.getNumberOfElements();
715 for (int ielem=0 ; ielem<nbElem; ielem++)
716 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
717 fres(ielem, jcomp)=_field->getArray()->getConstPointer()[jcomp+ielem*_field->getNumberOfComponents()]+f(ielem, jcomp);
721 //----------------------------------------------------------------------
723 Field::operator- ( const Field& f ) const
724 //----------------------------------------------------------------------
726 //if(f.getMesh().getMEDCouplingMesh() != _mesh.getMEDCouplingMesh())
727 //throw CdmathException("Field::operator- : Field subtraction requires identical meshes");
728 _mesh.getMEDCouplingMesh()->checkFastEquivalWith(f.getMesh().getMEDCouplingMesh(),1e-6);
729 if(f.getTypeOfField() != getTypeOfField())
730 throw CdmathException("Field::operator- : Field subtraction requires identical field types (CELLS, NODES or FACES");
731 if(f.getNumberOfComponents() != getNumberOfComponents())
732 throw CdmathException("Field::operator- : Field subtraction requires identical number of components");
734 Field fres(getName(),f.getTypeOfField(),f.getMesh(),f.getNumberOfComponents(),f.getTime());
735 int nbComp=f.getNumberOfComponents();
736 int nbElem=f.getNumberOfElements();
737 for (int ielem=0 ; ielem<nbElem; ielem++)
738 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
739 fres(ielem, jcomp)=_field->getArray()->getConstPointer()[jcomp+ielem*_field->getNumberOfComponents()]-f(ielem, jcomp);
743 //----------------------------------------------------------------------
745 Field::operator= ( const Field& f )
746 //----------------------------------------------------------------------
749 _typeField=f.getTypeOfField() ;
750 _numberOfComponents=f.getNumberOfComponents();
752 _fieldName=f.getName();
753 MCAuto<MEDCouplingFieldDouble> f1=f.getField()->deepCopy();
758 Field Field::deepCopy( ) const
760 Field F(getName(), getTypeOfField(), getMesh(), getNumberOfComponents(), getTime()) ;
761 MCAuto<MEDCouplingFieldDouble> f1=getField()->deepCopy();
762 F.setFieldByMEDCouplingFieldDouble(f1);
768 //----------------------------------------------------------------------
770 Field::operator+= ( const Field& f )
771 //----------------------------------------------------------------------
773 //if(f.getMesh().getMEDCouplingMesh() != _mesh.getMEDCouplingMesh())
774 //throw CdmathException("Field::operator+= : Field addition requires identical meshes");
775 _mesh.getMEDCouplingMesh()->checkFastEquivalWith(f.getMesh().getMEDCouplingMesh(),1e-6);
776 if(f.getTypeOfField() != getTypeOfField())
777 throw CdmathException("Field::operator+= : Field addition requires identical field types (CELLS, NODES or FACES");
778 if(f.getNumberOfComponents() != getNumberOfComponents())
779 throw CdmathException("Field::operator+= : Field addition requires identical number of components");
781 _field->setMesh(f.getField()->getMesh());
782 (*_field)+=(*f.getField());
786 //----------------------------------------------------------------------
788 Field::operator-= ( const Field& f )
789 //----------------------------------------------------------------------
791 //if(f.getMesh().getMEDCouplingMesh() != _mesh.getMEDCouplingMesh())
792 //throw CdmathException("Field::operator-= : Field subtraction requires identical meshes");
793 _mesh.getMEDCouplingMesh()->checkFastEquivalWith(f.getMesh().getMEDCouplingMesh(),1e-6);
794 if(f.getTypeOfField() != getTypeOfField())
795 throw CdmathException("Field::operator-= : Field subtraction requires identical field types (CELLS, NODES or FACES");
796 if(f.getNumberOfComponents() != getNumberOfComponents())
797 throw CdmathException("Field::operator-= : Field subtraction requires identical number of components");
799 _field->setMesh(f.getField()->getMesh());
800 (*_field)-=(*f.getField());
804 //----------------------------------------------------------------------
806 Field::operator*= ( double s )
807 //----------------------------------------------------------------------
809 int nbComp=getNumberOfComponents();
810 int nbElem=getNumberOfElements();
811 for (int i=0 ; i<nbComp ; i++)
812 for (int j=0 ; j<nbElem; j++)
813 _field->getArray()->getPointer()[i+j*_field->getNumberOfComponents()]*=s;
817 //----------------------------------------------------------------------
819 Field::operator/= ( double s )
820 //----------------------------------------------------------------------
822 int nbComp=getNumberOfComponents();
823 int nbElem=getNumberOfElements();
824 for (int i=0 ; i<nbComp ; i++)
825 for (int j=0 ; j<nbElem; j++)
826 _field->getArray()->getPointer()[i+j*_field->getNumberOfComponents()]/=s;
830 //----------------------------------------------------------------------
832 Field::operator-= ( double s )
833 //----------------------------------------------------------------------
835 int nbComp=getNumberOfComponents();
836 int nbElem=getNumberOfElements();
837 for (int i=0 ; i<nbComp ; i++)
838 for (int j=0 ; j<nbElem; j++)
839 _field->getArray()->getPointer()[i+j*_field->getNumberOfComponents()]-=s;
843 //----------------------------------------------------------------------
845 Field::operator+= ( double s )
846 //----------------------------------------------------------------------
848 int nbComp=getNumberOfComponents();
849 int nbElem=getNumberOfElements();
850 for (int i=0 ; i<nbComp ; i++)
851 for (int j=0 ; j<nbElem; j++)
852 _field->getArray()->getPointer()[i+j*_field->getNumberOfComponents()]+=s;
856 //----------------------------------------------------------------------
858 Field::writeVTK (std::string fileName, bool fromScratch) const
859 //----------------------------------------------------------------------
861 string fname=fileName+".pvd";
863 double time=_field->getTime(iter,order);
867 ofstream file(fname.c_str()) ;
868 file << "<VTKFile type=\"Collection\" version=\"0.1\" byte_order=\"LittleEndian\"><Collection>\n" ;
869 ostringstream numfile;
871 string filetmp=fileName+"_";
872 filetmp=filetmp+numfile.str();
873 string ret=_field->writeVTK(filetmp.c_str()) ;
874 file << "<DataSet timestep=\""<< time << "\" group=\"\" part=\"0\" file=\"" << ret << "\"/>\n" ;
875 file << "</Collection></VTKFile>\n" ;
880 ifstream file1(fname.c_str()) ;
882 getline(file1, contenus, '\0');
883 string to_remove="</Collection></VTKFile>";
884 size_t m = contenus.find(to_remove);
885 size_t n = contenus.find_first_of("\n", m + to_remove.length());
886 contenus.erase(m, n - m + 1);
888 ofstream file(fname.c_str()) ;
890 ostringstream numfile;
892 string filetmp=fileName+"_";
893 filetmp=filetmp+numfile.str();
894 string ret=_field->writeVTK(filetmp.c_str()) ;
895 file << "<DataSet timestep=\""<< time << "\" group=\"\" part=\"0\" file=\"" << ret << "\"/>\n" ;
896 file << "</Collection></VTKFile>\n" ;
901 //----------------------------------------------------------------------
903 Field::writeCSV ( const std::string fileName ) const
904 //----------------------------------------------------------------------
907 double time=_field->getTime(iter,order);
909 ostringstream numfile;
911 string filetmp=fileName+"_";
912 filetmp=filetmp+numfile.str();
913 filetmp=filetmp+".csv";
914 ofstream file(filetmp.c_str()) ;
915 int dim=_mesh.getSpaceDimension();
917 if (getTypeOfField()==CELLS)
918 nbElements=_mesh.getNumberOfCells();
920 nbElements=_mesh.getNumberOfNodes();
924 int nbCompo=getNumberOfComponents();
926 file << "x " << _field->getName() << endl;
930 for (int i=0;i<nbCompo;i++)
931 file << " " << _field->getName() << " Compo " << i+1 << " "<< getInfoOnComponent(i);
934 for (int i=0;i<nbElements;i++)
936 if (getTypeOfField()==CELLS)
937 file << _mesh.getCell(i).x() ;
939 file << _mesh.getNode(i).x() ;
940 for (int j=0;j<nbCompo;j++)
941 file << " " << getValues()[j+i*nbCompo] ;
946 int nbCompo=getNumberOfComponents();
948 file << "x y " << _field->getName() << endl;
952 for (int i=0;i<nbCompo;i++)
953 file << " " << _field->getName() << " Compo " << i+1<< " "<< getInfoOnComponent(i);
956 for (int i=0;i<nbElements;i++)
958 if (getTypeOfField()==CELLS)
959 file << _mesh.getCell(i).x() << " " << _mesh.getCell(i).y() ;
961 file << _mesh.getNode(i).x() << " " << _mesh.getNode(i).y() ;
962 for (int j=0;j<nbCompo;j++)
963 file << " " << getValues()[j+i*nbCompo] ;
968 int nbCompo=getNumberOfComponents();
970 file << "x y z " << _field->getName() << endl;
974 for (int i=0;i<nbCompo;i++)
975 file << " " << _field->getName() << " Compo " << i+1<< " "<< getInfoOnComponent(i);
978 for (int i=0;i<nbElements;i++)
980 if (getTypeOfField()==CELLS)
981 file << _mesh.getCell(i).x() << " " << _mesh.getCell(i).y() << " " << _mesh.getCell(i).z();
983 file << _mesh.getNode(i).x() << " " << _mesh.getNode(i).y() << " " << _mesh.getNode(i).z();
984 for (int j=0;j<nbCompo;j++)
985 file << " " << getValues()[j+i*nbCompo] ;
992 //----------------------------------------------------------------------
994 Field::writeMED ( const std::string fileName, bool fromScratch) const
995 //----------------------------------------------------------------------
997 string fname=fileName+".med";
999 MEDCoupling::WriteField(fname.c_str(),_field,fromScratch);
1001 MEDCoupling::WriteFieldUsingAlreadyWrittenMesh(fname.c_str(),_field);
1005 operator* (double value , const Field& field )
1007 Field fres(field.getName(),field.getTypeOfField(),field.getMesh(),field.getNumberOfComponents(),field.getTime());
1008 int nbComp=field.getNumberOfComponents();
1009 int nbElem=field.getNumberOfElements();
1010 for (int ielem=0 ; ielem<nbElem; ielem++)
1011 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
1012 fres(ielem, jcomp)=value*field(ielem, jcomp);
1017 operator* (const Field& field, double value )
1019 Field fres(field.getName(),field.getTypeOfField(),field.getMesh(),field.getNumberOfComponents(),field.getTime());
1020 int nbComp=field.getNumberOfComponents();
1021 int nbElem=field.getNumberOfElements();
1022 for (int ielem=0 ; ielem<nbElem; ielem++)
1023 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
1024 fres(ielem, jcomp)=value*field(ielem, jcomp);
1028 Field operator/ (const Field& field, double value)
1030 Field fres(field.getName(),field.getTypeOfField(),field.getMesh(),field.getNumberOfComponents(),field.getTime());
1031 int nbComp=field.getNumberOfComponents();
1032 int nbElem=field.getNumberOfElements();
1033 for (int ielem=0 ; ielem<nbElem; ielem++)
1034 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
1035 fres(ielem, jcomp)=field(ielem, jcomp)/value;
1040 Field::getValuesOnAllComponents(int elem) const
1042 Vector v(getNumberOfComponents());
1043 for(int i=0;i<getNumberOfComponents();i++)
1044 v[i]=(*this)(elem,i);
1049 Field::getValuesOnComponent(int compo) const
1051 Vector v(getNumberOfElements());
1052 for(int i=0;i<getNumberOfElements();i++)
1053 v[i]=(*this)(i,compo);
1057 std::ostream& operator<<(std::ostream& out, const Field& field )
1059 cout << "Field " << field.getName() << " : " << endl ;
1060 out<< field.getField().retn()->getArray()->repr();