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"
22 using namespace MEDCoupling;
26 //----------------------------------------------------------------------
27 Field::Field( EntityType typeField )
28 //----------------------------------------------------------------------
32 _numberOfComponents=0;
35 //----------------------------------------------------------------------
37 //----------------------------------------------------------------------
39 //std::cerr << "dtor Field, _field = " <<_field << std::endl;
40 //if (_field) _field->decrRef();
44 Field::Field(const std::string fieldName, EntityType type, const Mesh& mesh, int numberOfComponents, double time)
49 _numberOfComponents=numberOfComponents;
53 buildFieldMemoryStructure();
55 void Field::buildFieldMemoryStructure()
57 MEDCouplingUMesh* mu=_mesh.getMEDCouplingMesh()->buildUnstructured();
58 DataArrayDouble *array=DataArrayDouble::New();
59 if (_typeField==CELLS)
61 _field=MEDCouplingFieldDouble::New(ON_CELLS);
62 array->alloc(_mesh.getNumberOfCells(),_numberOfComponents);
64 }else if(_typeField==NODES)
66 _field=MEDCouplingFieldDouble::New(ON_NODES);
67 array->alloc(_mesh.getNumberOfNodes(),_numberOfComponents);
69 }else if(_typeField==FACES)
71 _field=MEDCouplingFieldDouble::New(ON_CELLS);
72 array->alloc(_mesh.getNumberOfFaces(),_numberOfComponents);
73 DataArrayIdType *desc=DataArrayIdType::New();
74 DataArrayIdType *descI=DataArrayIdType::New();
75 DataArrayIdType *revDesc=DataArrayIdType::New();
76 DataArrayIdType *revDescI=DataArrayIdType::New();
77 MEDCouplingUMesh *m3=mu->buildDescendingConnectivity(desc,descI,revDesc,revDescI);
85 throw CdmathException("Type of Field::Field() is not compatible");
87 _field->setName(_fieldName.c_str()) ;
88 _field->setArray(array);
89 _field->setTime(_time,0,0);
94 Field::Field( const std::string filename, EntityType type,
95 const std::string & fieldName,
96 int iteration, int order, int meshLevel,
97 int numberOfComponents, double time)
100 _mesh=Mesh(filename + ".med", meshLevel);
102 _numberOfComponents=numberOfComponents;
104 _fieldName=fieldName;
106 readFieldMed(filename, type, fieldName, iteration, order);
109 Field::Field(const std::string meshFileName, EntityType type, const std::vector<double> Vconstant,
110 const std::string & fieldName, int meshLevel, double time )
113 _mesh=Mesh(meshFileName + ".med", meshLevel);
115 _numberOfComponents=Vconstant.size();
117 _fieldName=fieldName;
119 buildFieldMemoryStructure();
121 int nbelem=_field->getNumberOfTuples();
122 int nbcomp=_field->getNumberOfComponents() ;
124 for (int ielem=0 ; ielem<nbelem; ielem++)
125 for (int jcomp=0 ; jcomp<nbcomp ; jcomp++)
126 _field->getArray()->getPointer()[jcomp+ielem*_field->getNumberOfComponents()]=Vconstant[jcomp];
128 Field::Field(const Mesh& M, EntityType type, const Vector Vconstant, const std::string & fieldName, double time)
133 _numberOfComponents=Vconstant.size();
135 _fieldName=fieldName;
137 buildFieldMemoryStructure();
139 int nbelem=_field->getNumberOfTuples();
140 int nbcomp=_field->getNumberOfComponents() ;
142 for (int ielem=0 ; ielem<nbelem; ielem++)
143 for (int jcomp=0 ; jcomp<nbcomp ; jcomp++)
144 _field->getArray()->getPointer()[jcomp+ielem*_field->getNumberOfComponents()]=Vconstant[jcomp];
146 Field::Field(const Mesh& M, EntityType type, const vector<double> Vconstant, const std::string & fieldName, double time)
151 _numberOfComponents=Vconstant.size();
153 _fieldName=fieldName;
155 buildFieldMemoryStructure();
157 int nbelem=_field->getNumberOfTuples();
158 int nbcomp=_field->getNumberOfComponents() ;
160 for (int ielem=0 ; ielem<nbelem; ielem++)
161 for (int jcomp=0 ; jcomp<nbcomp ; jcomp++)
162 _field->getArray()->getPointer()[jcomp+ielem*_field->getNumberOfComponents()]=Vconstant[jcomp];
164 Field::Field( int nDim, const vector<double> Vconstant, EntityType type,
165 double xmin, double xmax, int nx, string leftSide, string rightSide,
166 double ymin, double ymax, int ny, string backSide, string frontSide,
167 double zmin, double zmax, int nz, string bottomSide, string topSide,
168 const std::string & fieldName, double time,double epsilon)
172 _numberOfComponents=Vconstant.size();
174 _fieldName=fieldName;
178 _mesh=Mesh(xmin,xmax,nx);
181 _mesh=Mesh(xmin,xmax,nx,ymin,ymax,ny);
183 _mesh=Mesh(xmin,xmax,nx,ymin,ymax,ny,zmin,zmax,nz);
185 cout<<"Field::Field: Space dimension nDim should be between 1 and 3"<<endl;
186 throw CdmathException("Space dimension nDim should be between 1 and 3");
189 _mesh.setGroupAtPlan(xmax,0,epsilon,rightSide);
190 _mesh.setGroupAtPlan(xmin,0,epsilon,leftSide);
192 _mesh.setGroupAtPlan(ymax,1,epsilon,frontSide);
193 _mesh.setGroupAtPlan(ymin,1,epsilon,backSide);
196 _mesh.setGroupAtPlan(zmax,2,epsilon,topSide);
197 _mesh.setGroupAtPlan(zmin,2,epsilon,bottomSide);
201 buildFieldMemoryStructure();
203 int nbelem=_field->getNumberOfTuples();
204 int nbcomp=_field->getNumberOfComponents() ;
206 for (int ielem=0 ; ielem<nbelem; ielem++)
207 for (int jcomp=0 ; jcomp<nbcomp ; jcomp++)
208 _field->getArray()->getPointer()[jcomp+ielem*_field->getNumberOfComponents()]=Vconstant[jcomp];
210 Field::Field(const Mesh M, const Vector VV_Left, const Vector VV_Right, double disc_pos,
211 EntityType type, int direction, const std::string & fieldName, double time)
213 if (VV_Right.getNumberOfRows()!=VV_Left.getNumberOfRows())
214 throw CdmathException( "Field::Field: Vectors VV_Left and VV_Right have different sizes");
219 _numberOfComponents=VV_Left.getNumberOfRows();
221 _fieldName=fieldName;
224 buildFieldMemoryStructure();
226 int nbelem=_field->getNumberOfTuples();
227 int nbcomp=_field->getNumberOfComponents() ;
228 double component_value;
230 for (int j = 0; j < nbelem; j++) {
232 component_value=M.getCell(j).x();
233 else if(direction==1)
234 component_value=M.getCell(j).y();
235 else if(direction==2)
236 component_value=M.getCell(j).z();
238 throw CdmathException( "Field::Field: direction should be an integer between 0 and 2");
240 for (int i=0; i< nbcomp; i++)
241 if (component_value< disc_pos )
242 _field->getArray()->getPointer()[j+i*_field->getNumberOfComponents()] = VV_Left[i];
244 _field->getArray()->getPointer()[j+i*_field->getNumberOfComponents()] = VV_Right[i];
247 Field::Field( int nDim, const vector<double> VV_Left, vector<double> VV_Right,
248 double xstep, EntityType type,
249 double xmin, double xmax, int nx, string leftSide, string rightSide,
250 double ymin, double ymax, int ny, string backSide, string frontSide,
251 double zmin, double zmax, int nz, string bottomSide, string topSide,
252 int direction, const std::string & fieldName, double time, double epsilon)
254 if (VV_Right.size()!=VV_Left.size())
255 throw CdmathException( "Field::Field: Vectors VV_Left and VV_Right have different sizes");
258 M=Mesh(xmin,xmax,nx);
260 M=Mesh(xmin,xmax,nx,ymin,ymax,ny);
262 M=Mesh(xmin,xmax,nx,ymin,ymax,ny,zmin,zmax,nz);
264 throw CdmathException("Field::Field : Space dimension nDim should be between 1 and 3");
266 M.setGroupAtPlan(xmax,0,epsilon,rightSide);
267 M.setGroupAtPlan(xmin,0,epsilon,leftSide);
269 M.setGroupAtPlan(ymax,1,epsilon,frontSide);
270 M.setGroupAtPlan(ymin,1,epsilon,backSide);
273 M.setGroupAtPlan(zmax,2,epsilon,topSide);
274 M.setGroupAtPlan(zmin,2,epsilon,bottomSide);
276 Vector V_Left(VV_Left.size()), V_Right(VV_Right.size());
277 for(int i=0;i<VV_Left.size(); i++){
278 V_Left(i)=VV_Left[i];
279 V_Right(i)=VV_Right[i];
281 Field(M, V_Left, V_Right, xstep, type, direction, fieldName, time);
284 Field::Field(const Mesh M, const Vector Vin, const Vector Vout, double radius,
285 const Vector Center, EntityType type, const std::string & fieldName, double time)
287 if((Center.size()!=M.getSpaceDimension()) || (Vout.size() != Vin.size()) )
289 cout<< "Vout.size()= "<<Vout.size() << ", Vin.size()= "<<Vin.size()<<", Center.size()="<<Center.size()<<", M.getSpaceDim= "<< M.getSpaceDimension()<<endl;
290 throw CdmathException("Field::Field : Vector size error");
296 _numberOfComponents=Vout.size();
298 _fieldName=fieldName;
301 buildFieldMemoryStructure();
303 int nbelem=_field->getNumberOfTuples();
304 int nbcomp=_field->getNumberOfComponents() ;
306 int spaceDim=M.getSpaceDimension();
307 Vector currentPoint(spaceDim);
309 for(int i=0;i<nbelem;i++)
311 currentPoint(0)=M.getCell(i).x();
314 currentPoint(1)=M.getCell(i).y();
316 currentPoint(2)=M.getCell(i).z();
318 if((currentPoint-Center).norm()<radius)
319 for(int j=0;j<nbcomp;j++)
320 _field->getArray()->getPointer()[j+i*_field->getNumberOfComponents()]=Vin[j];
322 for(int j=0;j<nbcomp;j++)
323 _field->getArray()->getPointer()[j+i*_field->getNumberOfComponents()]=Vout[j];
327 MEDCoupling::DataArrayDouble * Field::getArray(){
328 return _field->getArray();
332 Field::readFieldMed( const std::string & fileNameRadical,
334 const std::string & fieldName,
339 * Reads the file fileNameRadical.med and creates a Field from it.
341 std::string completeFileName = fileNameRadical + ".med";
342 std::vector<std::string> fieldNames = MEDCoupling::GetAllFieldNames(completeFileName);
344 std::string attributedFieldName;
347 // Get the name of the right field that we will attribute to the Field.
348 if (fieldName == "") {
349 if (fieldNames.size() > 0)
350 attributedFieldName = fieldNames[0];
352 std::ostringstream message;
353 message << "No field in file " << completeFileName;
354 throw CdmathException(message.str().c_str());
358 for (; iField < fieldNames.size(); iField++)
359 if (fieldName == fieldNames[iField]) break;
361 if (iField < fieldNames.size())
362 attributedFieldName = fieldName;
364 std::ostringstream message;
365 message << "No field named " << fieldName << " in file " << completeFileName;
366 throw CdmathException(message.str().c_str());
370 // Get the name of the right mesh that we will attribute to the Field.
371 std::vector<std::string> meshNames
372 = MEDCoupling::GetMeshNamesOnField(completeFileName, attributedFieldName);
373 if (meshNames.size() == 0) {
374 std::ostringstream message;
375 message << "No mesh associated to " << fieldName
376 << " in file " << completeFileName;
377 throw CdmathException(message.str().c_str());
379 std::string attributedMeshName = meshNames[0];
382 MEDCoupling::TypeOfField medFieldType[3] = { ON_CELLS, ON_NODES, ON_CELLS };
385 _field = dynamic_cast< MEDCoupling::MEDCouplingFieldDouble * > (
386 MEDCoupling::ReadFieldCell( completeFileName,
387 attributedMeshName, 0,
388 attributedFieldName, iteration, order) );
391 _field = dynamic_cast< MEDCoupling::MEDCouplingFieldDouble * > (
392 MEDCoupling::ReadFieldNode( completeFileName,
393 attributedMeshName, 0,
394 attributedFieldName, iteration, order) );
397 _field = dynamic_cast< MEDCoupling::MEDCouplingFieldDouble * > (
398 MEDCoupling::ReadFieldCell( completeFileName,
399 attributedMeshName, -1,
400 attributedFieldName, iteration, order) );
407 Field::getNormEuclidean() const
409 Vector result(_numberOfComponents);
410 DoubleTab norm(_numberOfComponents,_field->magnitude()->getArray()->getConstPointer());
411 result.setValues(norm);
417 Field::max(int component) const
419 if( component >= getNumberOfComponents() || component < 0)
420 throw CdmathException("double Field::max() : component number should be between 0 and the field number of components");
422 double result=-1e100;
423 for(int i=0; i<getNumberOfElements() ; i++)
424 if( result < (*this)(i,component))
425 result = (*this)(i,component);
431 Field::min(int component) const
433 if( component >= getNumberOfComponents() || component < 0)
434 throw CdmathException("double Field::min() : component number should be between 0 and the field number of components");
437 for(int i=0; i<getNumberOfElements() ; i++)
438 if( result > (*this)(i,component))
439 result = (*this)(i,component);
445 Field::getInfoOnComponent(int icomp) const
447 return _field->getArray()->getInfoOnComponent(icomp);
451 Field::setInfoOnComponent(int icomp, string nameCompo)
453 _field.retn()->getArray()->setInfoOnComponent(icomp,nameCompo);
456 //----------------------------------------------------------------------
457 Field::Field( const Field & f )
458 //----------------------------------------------------------------------
461 MCAuto<MEDCouplingFieldDouble> f1=f.getField()->deepCopy();
463 _typeField=f.getTypeOfField();
466 //----------------------------------------------------------------------
467 MCAuto<MEDCouplingFieldDouble>
468 Field::getField ( void ) const
469 //----------------------------------------------------------------------
474 //----------------------------------------------------------------------
476 Field::setFieldByMEDCouplingFieldDouble ( const MEDCouplingFieldDouble* field )
477 //----------------------------------------------------------------------
479 MCAuto<MEDCouplingFieldDouble> ff=field->deepCopy();
483 //----------------------------------------------------------------------
485 Field::setFieldByDataArrayDouble ( const DataArrayDouble* array )
486 //----------------------------------------------------------------------
488 _field->setArray(const_cast<DataArrayDouble*>(array));
491 //----------------------------------------------------------------------
493 Field::integral ( ) const
494 //----------------------------------------------------------------------
496 int nbComp=_field->getNumberOfComponents();
497 double res[nbComp];//Pointer containing the inegral of each component
498 _field->integral(false,res);
499 Vector result(nbComp);//Vector containing the inegral of each component
501 for(int i=0; i<nbComp ; i++)
507 //----------------------------------------------------------------------
509 Field::integral ( int compId ) const
510 //----------------------------------------------------------------------
512 return _field->integral(compId, false);
515 //----------------------------------------------------------------------
517 Field::normL1 ( ) const
518 //----------------------------------------------------------------------
520 int nbComp=_field->getNumberOfComponents();
521 double res[nbComp];//Pointer containing the L1 norm of each component
523 Vector result(nbComp);//Vector containing the L1 norm of each component
525 for(int i=0; i<nbComp ; i++)
531 //----------------------------------------------------------------------
533 Field::normL2 ( ) const
534 //----------------------------------------------------------------------
536 int nbComp=_field->getNumberOfComponents();
537 double res[nbComp];//Pointer containing the L2 norm of each component
539 Vector result(nbComp);//Vector containing the L2 norm of each component
541 for(int i=0; i<nbComp ; i++)
547 //----------------------------------------------------------------------
549 Field::normMax ( ) const
550 //----------------------------------------------------------------------
552 int nbComp=_field->getNumberOfComponents();
553 double res[nbComp];//Pointer containing the L2 norm of each component
554 _field->normMax(res);
555 Vector result(nbComp);//Vector containing the L2 norm of each component
557 for(int i=0; i<nbComp ; i++)
564 Field::componentMax(Vector & Indices) const
566 int nbComp=_field->getNumberOfComponents();
567 int nbElems=getNumberOfElements();
569 Vector result(nbComp);//Vector containing the Linfinity norm of each component
571 for(int i=0; i<nbElems ; i++)
572 for(int j=0; j<nbComp ; j++)
573 if(fabs((*this)(i,j))>result(j))
575 result(j)=fabs((*this)(i,j));
581 //----------------------------------------------------------------------
583 Field::operator() ( int ielem )
584 //----------------------------------------------------------------------
586 if(ielem>_field->getNumberOfTuples() || ielem<0)
587 throw CdmathException("double& Field::operator(ielem) : ielem>number of values !");
588 return _field->getArray()->getPointer()[ielem*_field->getNumberOfComponents()];
591 //----------------------------------------------------------------------
593 Field::operator[] ( int ielem )
594 //----------------------------------------------------------------------
596 if(ielem>_field->getNumberOfTuples() || ielem<0)
597 throw CdmathException("double& Field::operator[ielem] : ielem>number of values !");
598 return _field->getArray()->getPointer()[ielem*_field->getNumberOfComponents()];
601 //----------------------------------------------------------------------
603 Field::operator() ( int ielem ) const
604 //----------------------------------------------------------------------
606 if(ielem>_field->getNumberOfTuples() || ielem<0)
607 throw CdmathException("double Field::operator(ielem) : ielem>number of values !");
608 return _field->getArray()->getConstPointer()[ielem*_field->getNumberOfComponents()];
611 //----------------------------------------------------------------------
613 Field::operator[] ( int ielem ) const
614 //----------------------------------------------------------------------
616 if(ielem>_field->getNumberOfTuples() || ielem<0)
617 throw CdmathException("double Field::operator[ielem] : ielem>number of values !");
618 return _field->getArray()->getConstPointer()[ielem*_field->getNumberOfComponents()];
621 //----------------------------------------------------------------------
623 Field::operator() ( int ielem, int jcomp )
624 //----------------------------------------------------------------------
626 if(ielem>_field->getNumberOfTuples() || jcomp>_field->getNumberOfComponents() || ielem<0 || jcomp<0)
627 throw CdmathException("double& Field::operator( int ielem, int jcomp ) : ielem>number of values or jcomp>number of components !");
628 return _field->getArray()->getPointer()[jcomp+ielem*_field->getNumberOfComponents()];
631 //----------------------------------------------------------------------
633 Field::operator() ( int ielem, int jcomp ) const
634 //----------------------------------------------------------------------
636 if(ielem>_field->getNumberOfTuples() || jcomp>_field->getNumberOfComponents() || ielem<0 || jcomp<0)
637 throw CdmathException("double Field::operator( int ielem, int jcomp ) : ielem>number of values or jcomp>number of components !");
638 return _field->getArray()->getConstPointer()[jcomp+ielem*_field->getNumberOfComponents()];
641 //----------------------------------------------------------------------
643 Field::setTime ( double time, int iter )
644 //----------------------------------------------------------------------
646 _field->setTime(time,iter,0.0);
648 //----------------------------------------------------------------------
650 Field::getTime ( void ) const
651 //----------------------------------------------------------------------
654 return _field->getTime(a,b);
657 //----------------------------------------------------------------------
659 Field::getNumberOfElements ( void ) const
660 //----------------------------------------------------------------------
662 return _field->getNumberOfTuples() ;
666 Field::getSpaceDimension( void ) const
668 return _mesh.getSpaceDimension() ;
671 //----------------------------------------------------------------------
673 Field::getNumberOfComponents ( void ) const
674 //----------------------------------------------------------------------
676 return _field->getNumberOfComponents() ;
679 //----------------------------------------------------------------------
681 Field::getValues ( void ) const
682 //----------------------------------------------------------------------
684 return _field->getArray()->getConstPointer() ;
687 //----------------------------------------------------------------------
689 Field::getValues ( Vector myVector ) const
690 //----------------------------------------------------------------------
692 if( myVector.size() != _field->getNumberOfTuples() * _field->getNumberOfComponents() )
693 throw CdmathException("Error : Field::getValues requires a vector having the same number of component as fiedl values");
695 const double * fieldValues = _field->getArray()->getConstPointer();
696 double * vectorValues = myVector.getValues().getValues();
698 memcpy(vectorValues, fieldValues,myVector.size()*sizeof(double)) ;
702 Field::setValues ( Vector myVector )
703 //----------------------------------------------------------------------
705 if( myVector.size() != _field->getNumberOfTuples() * _field->getNumberOfComponents() )
706 throw CdmathException("Error : Field::setValues requires a vector having the same number of component as fiedl values");
708 double * vectorValues = myVector.getValues().getValues();
710 setValues ( vectorValues, myVector.size() );
714 Field::setValues ( double * values, int nbValues )
716 double * fieldValues = _field->getArray()->getPointer() ;
717 memcpy(fieldValues,values,nbValues*sizeof(double)) ;
720 //----------------------------------------------------------------------
722 Field::getName ( void ) const
723 //----------------------------------------------------------------------
725 return _field->getName() ;
728 //----------------------------------------------------------------------
730 Field::getMesh ( void ) const
731 //----------------------------------------------------------------------
736 //----------------------------------------------------------------------
738 Field::getTypeOfField ( void ) const
739 //----------------------------------------------------------------------
744 //----------------------------------------------------------------------
746 Field::setName ( const string fieldName )
747 //----------------------------------------------------------------------
749 _field->setName(fieldName.c_str()) ;
753 //----------------------------------------------------------------------
755 Field::operator+ ( const Field& f ) const
756 //----------------------------------------------------------------------
758 //if(f.getMesh().getMEDCouplingMesh() != _mesh.getMEDCouplingMesh())
759 //throw CdmathException("Field::operator+ : Field addition requires identical meshes");
760 _mesh.getMEDCouplingMesh()->checkFastEquivalWith(f.getMesh().getMEDCouplingMesh(),1e-6);
761 if(f.getTypeOfField() != getTypeOfField())
762 throw CdmathException("Field::operator+ : Field addition requires identical field types (CELLS, NODES or FACES");
763 if(f.getNumberOfComponents() != getNumberOfComponents())
764 throw CdmathException("Field::operator+ : Field addition requires identical number of components");
766 Field fres(getName(),f.getTypeOfField(),f.getMesh(),f.getNumberOfComponents(),f.getTime());
767 int nbComp=f.getNumberOfComponents();
768 int nbElem=f.getNumberOfElements();
769 for (int ielem=0 ; ielem<nbElem; ielem++)
770 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
771 fres(ielem, jcomp)=_field->getArray()->getConstPointer()[jcomp+ielem*_field->getNumberOfComponents()]+f(ielem, jcomp);
775 //----------------------------------------------------------------------
777 Field::operator- ( const Field& f ) const
778 //----------------------------------------------------------------------
780 //if(f.getMesh().getMEDCouplingMesh() != _mesh.getMEDCouplingMesh())
781 //throw CdmathException("Field::operator- : Field subtraction requires identical meshes");
782 _mesh.getMEDCouplingMesh()->checkFastEquivalWith(f.getMesh().getMEDCouplingMesh(),1e-6);
783 if(f.getTypeOfField() != getTypeOfField())
784 throw CdmathException("Field::operator- : Field subtraction requires identical field types (CELLS, NODES or FACES");
785 if(f.getNumberOfComponents() != getNumberOfComponents())
786 throw CdmathException("Field::operator- : Field subtraction requires identical number of components");
788 Field fres(getName(),f.getTypeOfField(),f.getMesh(),f.getNumberOfComponents(),f.getTime());
789 int nbComp=f.getNumberOfComponents();
790 int nbElem=f.getNumberOfElements();
791 for (int ielem=0 ; ielem<nbElem; ielem++)
792 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
793 fres(ielem, jcomp)=_field->getArray()->getConstPointer()[jcomp+ielem*_field->getNumberOfComponents()]-f(ielem, jcomp);
797 //----------------------------------------------------------------------
799 Field::operator= ( const Field& f )
800 //----------------------------------------------------------------------
803 _typeField=f.getTypeOfField() ;
804 _numberOfComponents=f.getNumberOfComponents();
806 _fieldName=f.getName();
807 MCAuto<MEDCouplingFieldDouble> f1=f.getField()->deepCopy();
812 Field Field::deepCopy( ) const
814 Field F(getName(), getTypeOfField(), getMesh(), getNumberOfComponents(), getTime()) ;
815 MCAuto<MEDCouplingFieldDouble> f1=getField()->deepCopy();
816 F.setFieldByMEDCouplingFieldDouble(f1);
822 //----------------------------------------------------------------------
824 Field::operator+= ( const Field& f )
825 //----------------------------------------------------------------------
827 //if(f.getMesh().getMEDCouplingMesh() != _mesh.getMEDCouplingMesh())
828 //throw CdmathException("Field::operator+= : Field addition requires identical meshes");
829 _mesh.getMEDCouplingMesh()->checkFastEquivalWith(f.getMesh().getMEDCouplingMesh(),1e-6);
830 if(f.getTypeOfField() != getTypeOfField())
831 throw CdmathException("Field::operator+= : Field addition requires identical field types (CELLS, NODES or FACES");
832 if(f.getNumberOfComponents() != getNumberOfComponents())
833 throw CdmathException("Field::operator+= : Field addition requires identical number of components");
835 _field->setMesh(f.getField()->getMesh());
836 (*_field)+=(*f.getField());
840 //----------------------------------------------------------------------
842 Field::operator-= ( const Field& f )
843 //----------------------------------------------------------------------
845 //if(f.getMesh().getMEDCouplingMesh() != _mesh.getMEDCouplingMesh())
846 //throw CdmathException("Field::operator-= : Field subtraction requires identical meshes");
847 _mesh.getMEDCouplingMesh()->checkFastEquivalWith(f.getMesh().getMEDCouplingMesh(),1e-6);
848 if(f.getTypeOfField() != getTypeOfField())
849 throw CdmathException("Field::operator-= : Field subtraction requires identical field types (CELLS, NODES or FACES");
850 if(f.getNumberOfComponents() != getNumberOfComponents())
851 throw CdmathException("Field::operator-= : Field subtraction requires identical number of components");
853 _field->setMesh(f.getField()->getMesh());
854 (*_field)-=(*f.getField());
858 //----------------------------------------------------------------------
860 Field::operator*= ( double s )
861 //----------------------------------------------------------------------
863 int nbComp=getNumberOfComponents();
864 int nbElem=getNumberOfElements();
865 for (int i=0 ; i<nbComp ; i++)
866 for (int j=0 ; j<nbElem; j++)
867 _field->getArray()->getPointer()[i+j*_field->getNumberOfComponents()]*=s;
871 //----------------------------------------------------------------------
873 Field::operator/= ( double s )
874 //----------------------------------------------------------------------
876 int nbComp=getNumberOfComponents();
877 int nbElem=getNumberOfElements();
878 for (int i=0 ; i<nbComp ; i++)
879 for (int j=0 ; j<nbElem; j++)
880 _field->getArray()->getPointer()[i+j*_field->getNumberOfComponents()]/=s;
884 //----------------------------------------------------------------------
886 Field::operator-= ( double s )
887 //----------------------------------------------------------------------
889 int nbComp=getNumberOfComponents();
890 int nbElem=getNumberOfElements();
891 for (int i=0 ; i<nbComp ; i++)
892 for (int j=0 ; j<nbElem; j++)
893 _field->getArray()->getPointer()[i+j*_field->getNumberOfComponents()]-=s;
897 //----------------------------------------------------------------------
899 Field::operator+= ( double s )
900 //----------------------------------------------------------------------
902 int nbComp=getNumberOfComponents();
903 int nbElem=getNumberOfElements();
904 for (int i=0 ; i<nbComp ; i++)
905 for (int j=0 ; j<nbElem; j++)
906 _field->getArray()->getPointer()[i+j*_field->getNumberOfComponents()]+=s;
910 //----------------------------------------------------------------------
912 Field::writeVTK (std::string fileName, bool fromScratch) const
913 //----------------------------------------------------------------------
915 string fname=fileName+".pvd";
917 double time=_field->getTime(iter,order);
921 ofstream file(fname.c_str()) ;
922 file << "<VTKFile type=\"Collection\" version=\"0.1\" byte_order=\"LittleEndian\"><Collection>\n" ;
923 ostringstream numfile;
925 string filetmp=fileName+"_";
926 filetmp=filetmp+numfile.str();
927 string ret=_field->writeVTK(filetmp.c_str()) ;
928 file << "<DataSet timestep=\""<< time << "\" group=\"\" part=\"0\" file=\"" << ret << "\"/>\n" ;
929 file << "</Collection></VTKFile>\n" ;
934 ifstream file1(fname.c_str()) ;
936 getline(file1, contenus, '\0');
937 string to_remove="</Collection></VTKFile>";
938 size_t m = contenus.find(to_remove);
939 size_t n = contenus.find_first_of("\n", m + to_remove.length());
940 contenus.erase(m, n - m + 1);
942 ofstream file(fname.c_str()) ;
944 ostringstream numfile;
946 string filetmp=fileName+"_";
947 filetmp=filetmp+numfile.str();
948 string ret=_field->writeVTK(filetmp.c_str()) ;
949 file << "<DataSet timestep=\""<< time << "\" group=\"\" part=\"0\" file=\"" << ret << "\"/>\n" ;
950 file << "</Collection></VTKFile>\n" ;
955 //----------------------------------------------------------------------
957 Field::writeCSV ( const std::string fileName ) const
958 //----------------------------------------------------------------------
961 double time=_field->getTime(iter,order);
963 ostringstream numfile;
965 string filetmp=fileName+"_";
966 filetmp=filetmp+numfile.str();
967 filetmp=filetmp+".csv";
968 ofstream file(filetmp.c_str()) ;
969 int dim=_mesh.getSpaceDimension();
971 if (getTypeOfField()==CELLS)
972 nbElements=_mesh.getNumberOfCells();
974 nbElements=_mesh.getNumberOfNodes();
978 int nbCompo=getNumberOfComponents();
980 file << "x " << _field->getName() << endl;
984 for (int i=0;i<nbCompo;i++)
985 file << " " << _field->getName() << " Compo " << i+1 << " "<< getInfoOnComponent(i);
988 for (int i=0;i<nbElements;i++)
990 if (getTypeOfField()==CELLS)
991 file << _mesh.getCell(i).x() ;
993 file << _mesh.getNode(i).x() ;
994 for (int j=0;j<nbCompo;j++)
995 file << " " << getValues()[j+i*nbCompo] ;
1000 int nbCompo=getNumberOfComponents();
1002 file << "x y " << _field->getName() << endl;
1006 for (int i=0;i<nbCompo;i++)
1007 file << " " << _field->getName() << " Compo " << i+1<< " "<< getInfoOnComponent(i);
1010 for (int i=0;i<nbElements;i++)
1012 if (getTypeOfField()==CELLS)
1013 file << _mesh.getCell(i).x() << " " << _mesh.getCell(i).y() ;
1015 file << _mesh.getNode(i).x() << " " << _mesh.getNode(i).y() ;
1016 for (int j=0;j<nbCompo;j++)
1017 file << " " << getValues()[j+i*nbCompo] ;
1022 int nbCompo=getNumberOfComponents();
1024 file << "x y z " << _field->getName() << endl;
1028 for (int i=0;i<nbCompo;i++)
1029 file << " " << _field->getName() << " Compo " << i+1<< " "<< getInfoOnComponent(i);
1032 for (int i=0;i<nbElements;i++)
1034 if (getTypeOfField()==CELLS)
1035 file << _mesh.getCell(i).x() << " " << _mesh.getCell(i).y() << " " << _mesh.getCell(i).z();
1037 file << _mesh.getNode(i).x() << " " << _mesh.getNode(i).y() << " " << _mesh.getNode(i).z();
1038 for (int j=0;j<nbCompo;j++)
1039 file << " " << getValues()[j+i*nbCompo] ;
1046 //----------------------------------------------------------------------
1048 Field::writeMED ( const std::string fileName, bool fromScratch) const
1049 //----------------------------------------------------------------------
1051 string fname=fileName+".med";
1053 MEDCoupling::WriteField(fname.c_str(),_field,fromScratch);
1055 MEDCoupling::WriteFieldUsingAlreadyWrittenMesh(fname.c_str(),_field);
1059 operator* (double value , const Field& field )
1061 Field fres(field.getName(),field.getTypeOfField(),field.getMesh(),field.getNumberOfComponents(),field.getTime());
1062 int nbComp=field.getNumberOfComponents();
1063 int nbElem=field.getNumberOfElements();
1064 for (int ielem=0 ; ielem<nbElem; ielem++)
1065 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
1066 fres(ielem, jcomp)=value*field(ielem, jcomp);
1071 operator* (const Field& field, double value )
1073 Field fres(field.getName(),field.getTypeOfField(),field.getMesh(),field.getNumberOfComponents(),field.getTime());
1074 int nbComp=field.getNumberOfComponents();
1075 int nbElem=field.getNumberOfElements();
1076 for (int ielem=0 ; ielem<nbElem; ielem++)
1077 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
1078 fres(ielem, jcomp)=value*field(ielem, jcomp);
1082 Field operator/ (const Field& field, double value)
1084 Field fres(field.getName(),field.getTypeOfField(),field.getMesh(),field.getNumberOfComponents(),field.getTime());
1085 int nbComp=field.getNumberOfComponents();
1086 int nbElem=field.getNumberOfElements();
1087 for (int ielem=0 ; ielem<nbElem; ielem++)
1088 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
1089 fres(ielem, jcomp)=field(ielem, jcomp)/value;
1094 Field::getValuesOnAllComponents(int elem) const
1096 Vector v(getNumberOfComponents());
1097 for(int i=0;i<getNumberOfComponents();i++)
1098 v[i]=(*this)(elem,i);
1103 Field::getValuesOnComponent(int compo) const
1105 Vector v(getNumberOfElements());
1106 for(int i=0;i<getNumberOfElements();i++)
1107 v[i]=(*this)(i,compo);
1111 std::vector< double >
1112 Field::getFieldValues(int compo) const
1114 std::vector< double > v(getNumberOfElements());
1115 for(int i=0;i<getNumberOfElements();i++)
1116 v[i]=(*this)(i,compo);
1120 std::ostream& operator<<(std::ostream& out, const Field& field )
1122 cout << "Field " << field.getName() << " : " << endl ;
1123 out<< field.getField().retn()->getArray()->repr();