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 Vector result(_numberOfComponents);
408 DoubleTab norm(_numberOfComponents,_field->magnitude()->getArray()->getConstPointer());
409 result.setValues(norm);
415 Field::max(int component) const
417 if( component >= getNumberOfComponents() || component < 0)
418 throw CdmathException("double Field::max() : component number should be between 0 and the field number of components");
420 double result=-1e100;
421 for(int i=0; i<getNumberOfElements() ; i++)
422 if( result < (*this)(i,component))
423 result = (*this)(i,component);
429 Field::min(int component) const
431 if( component >= getNumberOfComponents() || component < 0)
432 throw CdmathException("double Field::min() : component number should be between 0 and the field number of components");
435 for(int i=0; i<getNumberOfElements() ; i++)
436 if( result > (*this)(i,component))
437 result = (*this)(i,component);
443 Field::getInfoOnComponent(int icomp) const
445 return _field->getArray()->getInfoOnComponent(icomp);
449 Field::setInfoOnComponent(int icomp, string nameCompo)
451 _field.retn()->getArray()->setInfoOnComponent(icomp,nameCompo);
454 //----------------------------------------------------------------------
455 Field::Field( const Field & f )
456 //----------------------------------------------------------------------
459 MCAuto<MEDCouplingFieldDouble> f1=f.getField()->deepCopy();
461 _typeField=f.getTypeOfField();
464 //----------------------------------------------------------------------
465 MCAuto<MEDCouplingFieldDouble>
466 Field::getField ( void ) const
467 //----------------------------------------------------------------------
472 //----------------------------------------------------------------------
474 Field::setFieldByMEDCouplingFieldDouble ( const MEDCouplingFieldDouble* field )
475 //----------------------------------------------------------------------
477 MCAuto<MEDCouplingFieldDouble> ff=field->deepCopy();
481 //----------------------------------------------------------------------
483 Field::setFieldByDataArrayDouble ( const DataArrayDouble* array )
484 //----------------------------------------------------------------------
486 _field->setArray(const_cast<DataArrayDouble*>(array));
489 //----------------------------------------------------------------------
491 Field::integral ( ) const
492 //----------------------------------------------------------------------
494 int nbComp=_field->getNumberOfComponents();
495 double res[nbComp];//Pointer containing the inegral of each component
496 _field->integral(false,res);
497 Vector result(nbComp);//Vector containing the inegral of each component
499 for(int i=0; i<nbComp ; i++)
505 //----------------------------------------------------------------------
507 Field::integral ( int compId ) const
508 //----------------------------------------------------------------------
510 return _field->integral(compId, false);
513 //----------------------------------------------------------------------
515 Field::normL1 ( ) const
516 //----------------------------------------------------------------------
518 int nbComp=_field->getNumberOfComponents();
519 double res[nbComp];//Pointer containing the L1 norm of each component
521 Vector result(nbComp);//Vector containing the L1 norm of each component
523 for(int i=0; i<nbComp ; i++)
529 //----------------------------------------------------------------------
531 Field::normL2 ( ) const
532 //----------------------------------------------------------------------
534 int nbComp=_field->getNumberOfComponents();
535 double res[nbComp];//Pointer containing the L2 norm of each component
537 Vector result(nbComp);//Vector containing the L2 norm of each component
539 for(int i=0; i<nbComp ; i++)
545 //----------------------------------------------------------------------
547 Field::normMax ( ) const
548 //----------------------------------------------------------------------
550 int nbComp=_field->getNumberOfComponents();
551 double res[nbComp];//Pointer containing the L2 norm of each component
552 _field->normMax(res);
553 Vector result(nbComp);//Vector containing the L2 norm of each component
555 for(int i=0; i<nbComp ; i++)
562 Field::componentMax(Vector & Indices) const
564 int nbComp=_field->getNumberOfComponents();
565 int nbElems=getNumberOfElements();
567 Vector result(nbComp);//Vector containing the Linfinity norm of each component
569 for(int i=0; i<nbElems ; i++)
570 for(int j=0; j<nbComp ; j++)
571 if(fabs((*this)(i,j))>result(j))
573 result(j)=fabs((*this)(i,j));
579 //----------------------------------------------------------------------
581 Field::operator() ( int ielem )
582 //----------------------------------------------------------------------
584 if(ielem>_field->getNumberOfTuples() || ielem<0)
585 throw CdmathException("double& Field::operator(ielem) : ielem>number of values !");
586 return _field->getArray()->getPointer()[ielem*_field->getNumberOfComponents()];
589 //----------------------------------------------------------------------
591 Field::operator[] ( int ielem )
592 //----------------------------------------------------------------------
594 if(ielem>_field->getNumberOfTuples() || ielem<0)
595 throw CdmathException("double& Field::operator[ielem] : ielem>number of values !");
596 return _field->getArray()->getPointer()[ielem*_field->getNumberOfComponents()];
599 //----------------------------------------------------------------------
601 Field::operator() ( int ielem ) const
602 //----------------------------------------------------------------------
604 if(ielem>_field->getNumberOfTuples() || ielem<0)
605 throw CdmathException("double Field::operator(ielem) : ielem>number of values !");
606 return _field->getArray()->getConstPointer()[ielem*_field->getNumberOfComponents()];
609 //----------------------------------------------------------------------
611 Field::operator[] ( int ielem ) const
612 //----------------------------------------------------------------------
614 if(ielem>_field->getNumberOfTuples() || ielem<0)
615 throw CdmathException("double Field::operator[ielem] : ielem>number of values !");
616 return _field->getArray()->getConstPointer()[ielem*_field->getNumberOfComponents()];
619 //----------------------------------------------------------------------
621 Field::operator() ( int ielem, int jcomp )
622 //----------------------------------------------------------------------
624 if(ielem>_field->getNumberOfTuples() || jcomp>_field->getNumberOfComponents() || ielem<0 || jcomp<0)
625 throw CdmathException("double& Field::operator( int ielem, int jcomp ) : ielem>number of values or jcomp>number of components !");
626 return _field->getArray()->getPointer()[jcomp+ielem*_field->getNumberOfComponents()];
629 //----------------------------------------------------------------------
631 Field::operator() ( int ielem, int jcomp ) const
632 //----------------------------------------------------------------------
634 if(ielem>_field->getNumberOfTuples() || jcomp>_field->getNumberOfComponents() || ielem<0 || jcomp<0)
635 throw CdmathException("double Field::operator( int ielem, int jcomp ) : ielem>number of values or jcomp>number of components !");
636 return _field->getArray()->getConstPointer()[jcomp+ielem*_field->getNumberOfComponents()];
639 //----------------------------------------------------------------------
641 Field::setTime ( double time, int iter )
642 //----------------------------------------------------------------------
644 _field->setTime(time,iter,0.0);
646 //----------------------------------------------------------------------
648 Field::getTime ( void ) const
649 //----------------------------------------------------------------------
652 return _field->getTime(a,b);
655 //----------------------------------------------------------------------
657 Field::getNumberOfElements ( void ) const
658 //----------------------------------------------------------------------
660 return _field->getNumberOfTuples() ;
664 Field::getSpaceDimension( void ) const
666 return _mesh.getSpaceDimension() ;
669 //----------------------------------------------------------------------
671 Field::getNumberOfComponents ( void ) const
672 //----------------------------------------------------------------------
674 return _field->getNumberOfComponents() ;
677 //----------------------------------------------------------------------
679 Field::getValues ( void ) const
680 //----------------------------------------------------------------------
682 return _field->getArray()->getConstPointer() ;
685 //----------------------------------------------------------------------
687 Field::getName ( void ) const
688 //----------------------------------------------------------------------
690 return _field->getName() ;
693 //----------------------------------------------------------------------
695 Field::getMesh ( void ) const
696 //----------------------------------------------------------------------
701 //----------------------------------------------------------------------
703 Field::getTypeOfField ( void ) const
704 //----------------------------------------------------------------------
709 //----------------------------------------------------------------------
711 Field::setName ( const string fieldName )
712 //----------------------------------------------------------------------
714 _field->setName(fieldName.c_str()) ;
718 //----------------------------------------------------------------------
720 Field::operator+ ( const Field& f ) const
721 //----------------------------------------------------------------------
723 //if(f.getMesh().getMEDCouplingMesh() != _mesh.getMEDCouplingMesh())
724 //throw CdmathException("Field::operator+ : Field addition requires identical meshes");
725 _mesh.getMEDCouplingMesh()->checkFastEquivalWith(f.getMesh().getMEDCouplingMesh(),1e-6);
726 if(f.getTypeOfField() != getTypeOfField())
727 throw CdmathException("Field::operator+ : Field addition requires identical field types (CELLS, NODES or FACES");
728 if(f.getNumberOfComponents() != getNumberOfComponents())
729 throw CdmathException("Field::operator+ : Field addition requires identical number of components");
731 Field fres(getName(),f.getTypeOfField(),f.getMesh(),f.getNumberOfComponents(),f.getTime());
732 int nbComp=f.getNumberOfComponents();
733 int nbElem=f.getNumberOfElements();
734 for (int ielem=0 ; ielem<nbElem; ielem++)
735 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
736 fres(ielem, jcomp)=_field->getArray()->getConstPointer()[jcomp+ielem*_field->getNumberOfComponents()]+f(ielem, jcomp);
740 //----------------------------------------------------------------------
742 Field::operator- ( const Field& f ) const
743 //----------------------------------------------------------------------
745 //if(f.getMesh().getMEDCouplingMesh() != _mesh.getMEDCouplingMesh())
746 //throw CdmathException("Field::operator- : Field subtraction requires identical meshes");
747 _mesh.getMEDCouplingMesh()->checkFastEquivalWith(f.getMesh().getMEDCouplingMesh(),1e-6);
748 if(f.getTypeOfField() != getTypeOfField())
749 throw CdmathException("Field::operator- : Field subtraction requires identical field types (CELLS, NODES or FACES");
750 if(f.getNumberOfComponents() != getNumberOfComponents())
751 throw CdmathException("Field::operator- : Field subtraction requires identical number of components");
753 Field fres(getName(),f.getTypeOfField(),f.getMesh(),f.getNumberOfComponents(),f.getTime());
754 int nbComp=f.getNumberOfComponents();
755 int nbElem=f.getNumberOfElements();
756 for (int ielem=0 ; ielem<nbElem; ielem++)
757 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
758 fres(ielem, jcomp)=_field->getArray()->getConstPointer()[jcomp+ielem*_field->getNumberOfComponents()]-f(ielem, jcomp);
762 //----------------------------------------------------------------------
764 Field::operator= ( const Field& f )
765 //----------------------------------------------------------------------
768 _typeField=f.getTypeOfField() ;
769 _numberOfComponents=f.getNumberOfComponents();
771 _fieldName=f.getName();
772 MCAuto<MEDCouplingFieldDouble> f1=f.getField()->deepCopy();
777 Field Field::deepCopy( ) const
779 Field F(getName(), getTypeOfField(), getMesh(), getNumberOfComponents(), getTime()) ;
780 MCAuto<MEDCouplingFieldDouble> f1=getField()->deepCopy();
781 F.setFieldByMEDCouplingFieldDouble(f1);
787 //----------------------------------------------------------------------
789 Field::operator+= ( const Field& f )
790 //----------------------------------------------------------------------
792 //if(f.getMesh().getMEDCouplingMesh() != _mesh.getMEDCouplingMesh())
793 //throw CdmathException("Field::operator+= : Field addition requires identical meshes");
794 _mesh.getMEDCouplingMesh()->checkFastEquivalWith(f.getMesh().getMEDCouplingMesh(),1e-6);
795 if(f.getTypeOfField() != getTypeOfField())
796 throw CdmathException("Field::operator+= : Field addition requires identical field types (CELLS, NODES or FACES");
797 if(f.getNumberOfComponents() != getNumberOfComponents())
798 throw CdmathException("Field::operator+= : Field addition requires identical number of components");
800 _field->setMesh(f.getField()->getMesh());
801 (*_field)+=(*f.getField());
805 //----------------------------------------------------------------------
807 Field::operator-= ( const Field& f )
808 //----------------------------------------------------------------------
810 //if(f.getMesh().getMEDCouplingMesh() != _mesh.getMEDCouplingMesh())
811 //throw CdmathException("Field::operator-= : Field subtraction requires identical meshes");
812 _mesh.getMEDCouplingMesh()->checkFastEquivalWith(f.getMesh().getMEDCouplingMesh(),1e-6);
813 if(f.getTypeOfField() != getTypeOfField())
814 throw CdmathException("Field::operator-= : Field subtraction requires identical field types (CELLS, NODES or FACES");
815 if(f.getNumberOfComponents() != getNumberOfComponents())
816 throw CdmathException("Field::operator-= : Field subtraction requires identical number of components");
818 _field->setMesh(f.getField()->getMesh());
819 (*_field)-=(*f.getField());
823 //----------------------------------------------------------------------
825 Field::operator*= ( double s )
826 //----------------------------------------------------------------------
828 int nbComp=getNumberOfComponents();
829 int nbElem=getNumberOfElements();
830 for (int i=0 ; i<nbComp ; i++)
831 for (int j=0 ; j<nbElem; j++)
832 _field->getArray()->getPointer()[i+j*_field->getNumberOfComponents()]*=s;
836 //----------------------------------------------------------------------
838 Field::operator/= ( double s )
839 //----------------------------------------------------------------------
841 int nbComp=getNumberOfComponents();
842 int nbElem=getNumberOfElements();
843 for (int i=0 ; i<nbComp ; i++)
844 for (int j=0 ; j<nbElem; j++)
845 _field->getArray()->getPointer()[i+j*_field->getNumberOfComponents()]/=s;
849 //----------------------------------------------------------------------
851 Field::operator-= ( double s )
852 //----------------------------------------------------------------------
854 int nbComp=getNumberOfComponents();
855 int nbElem=getNumberOfElements();
856 for (int i=0 ; i<nbComp ; i++)
857 for (int j=0 ; j<nbElem; j++)
858 _field->getArray()->getPointer()[i+j*_field->getNumberOfComponents()]-=s;
862 //----------------------------------------------------------------------
864 Field::operator+= ( double s )
865 //----------------------------------------------------------------------
867 int nbComp=getNumberOfComponents();
868 int nbElem=getNumberOfElements();
869 for (int i=0 ; i<nbComp ; i++)
870 for (int j=0 ; j<nbElem; j++)
871 _field->getArray()->getPointer()[i+j*_field->getNumberOfComponents()]+=s;
875 //----------------------------------------------------------------------
877 Field::writeVTK (std::string fileName, bool fromScratch) const
878 //----------------------------------------------------------------------
880 string fname=fileName+".pvd";
882 double time=_field->getTime(iter,order);
886 ofstream file(fname.c_str()) ;
887 file << "<VTKFile type=\"Collection\" version=\"0.1\" byte_order=\"LittleEndian\"><Collection>\n" ;
888 ostringstream numfile;
890 string filetmp=fileName+"_";
891 filetmp=filetmp+numfile.str();
892 string ret=_field->writeVTK(filetmp.c_str()) ;
893 file << "<DataSet timestep=\""<< time << "\" group=\"\" part=\"0\" file=\"" << ret << "\"/>\n" ;
894 file << "</Collection></VTKFile>\n" ;
899 ifstream file1(fname.c_str()) ;
901 getline(file1, contenus, '\0');
902 string to_remove="</Collection></VTKFile>";
903 size_t m = contenus.find(to_remove);
904 size_t n = contenus.find_first_of("\n", m + to_remove.length());
905 contenus.erase(m, n - m + 1);
907 ofstream file(fname.c_str()) ;
909 ostringstream numfile;
911 string filetmp=fileName+"_";
912 filetmp=filetmp+numfile.str();
913 string ret=_field->writeVTK(filetmp.c_str()) ;
914 file << "<DataSet timestep=\""<< time << "\" group=\"\" part=\"0\" file=\"" << ret << "\"/>\n" ;
915 file << "</Collection></VTKFile>\n" ;
920 //----------------------------------------------------------------------
922 Field::writeCSV ( const std::string fileName ) const
923 //----------------------------------------------------------------------
926 double time=_field->getTime(iter,order);
928 ostringstream numfile;
930 string filetmp=fileName+"_";
931 filetmp=filetmp+numfile.str();
932 filetmp=filetmp+".csv";
933 ofstream file(filetmp.c_str()) ;
934 int dim=_mesh.getSpaceDimension();
936 if (getTypeOfField()==CELLS)
937 nbElements=_mesh.getNumberOfCells();
939 nbElements=_mesh.getNumberOfNodes();
943 int nbCompo=getNumberOfComponents();
945 file << "x " << _field->getName() << endl;
949 for (int i=0;i<nbCompo;i++)
950 file << " " << _field->getName() << " Compo " << i+1 << " "<< getInfoOnComponent(i);
953 for (int i=0;i<nbElements;i++)
955 if (getTypeOfField()==CELLS)
956 file << _mesh.getCell(i).x() ;
958 file << _mesh.getNode(i).x() ;
959 for (int j=0;j<nbCompo;j++)
960 file << " " << getValues()[j+i*nbCompo] ;
965 int nbCompo=getNumberOfComponents();
967 file << "x y " << _field->getName() << endl;
971 for (int i=0;i<nbCompo;i++)
972 file << " " << _field->getName() << " Compo " << i+1<< " "<< getInfoOnComponent(i);
975 for (int i=0;i<nbElements;i++)
977 if (getTypeOfField()==CELLS)
978 file << _mesh.getCell(i).x() << " " << _mesh.getCell(i).y() ;
980 file << _mesh.getNode(i).x() << " " << _mesh.getNode(i).y() ;
981 for (int j=0;j<nbCompo;j++)
982 file << " " << getValues()[j+i*nbCompo] ;
987 int nbCompo=getNumberOfComponents();
989 file << "x y z " << _field->getName() << endl;
993 for (int i=0;i<nbCompo;i++)
994 file << " " << _field->getName() << " Compo " << i+1<< " "<< getInfoOnComponent(i);
997 for (int i=0;i<nbElements;i++)
999 if (getTypeOfField()==CELLS)
1000 file << _mesh.getCell(i).x() << " " << _mesh.getCell(i).y() << " " << _mesh.getCell(i).z();
1002 file << _mesh.getNode(i).x() << " " << _mesh.getNode(i).y() << " " << _mesh.getNode(i).z();
1003 for (int j=0;j<nbCompo;j++)
1004 file << " " << getValues()[j+i*nbCompo] ;
1011 //----------------------------------------------------------------------
1013 Field::writeMED ( const std::string fileName, bool fromScratch) const
1014 //----------------------------------------------------------------------
1016 string fname=fileName+".med";
1018 MEDCoupling::WriteField(fname.c_str(),_field,fromScratch);
1020 MEDCoupling::WriteFieldUsingAlreadyWrittenMesh(fname.c_str(),_field);
1024 operator* (double value , const Field& field )
1026 Field fres(field.getName(),field.getTypeOfField(),field.getMesh(),field.getNumberOfComponents(),field.getTime());
1027 int nbComp=field.getNumberOfComponents();
1028 int nbElem=field.getNumberOfElements();
1029 for (int ielem=0 ; ielem<nbElem; ielem++)
1030 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
1031 fres(ielem, jcomp)=value*field(ielem, jcomp);
1036 operator* (const Field& field, double value )
1038 Field fres(field.getName(),field.getTypeOfField(),field.getMesh(),field.getNumberOfComponents(),field.getTime());
1039 int nbComp=field.getNumberOfComponents();
1040 int nbElem=field.getNumberOfElements();
1041 for (int ielem=0 ; ielem<nbElem; ielem++)
1042 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
1043 fres(ielem, jcomp)=value*field(ielem, jcomp);
1047 Field operator/ (const Field& field, double value)
1049 Field fres(field.getName(),field.getTypeOfField(),field.getMesh(),field.getNumberOfComponents(),field.getTime());
1050 int nbComp=field.getNumberOfComponents();
1051 int nbElem=field.getNumberOfElements();
1052 for (int ielem=0 ; ielem<nbElem; ielem++)
1053 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
1054 fres(ielem, jcomp)=field(ielem, jcomp)/value;
1059 Field::getValuesOnAllComponents(int elem) const
1061 Vector v(getNumberOfComponents());
1062 for(int i=0;i<getNumberOfComponents();i++)
1063 v[i]=(*this)(elem,i);
1068 Field::getValuesOnComponent(int compo) const
1070 Vector v(getNumberOfElements());
1071 for(int i=0;i<getNumberOfElements();i++)
1072 v[i]=(*this)(i,compo);
1076 std::ostream& operator<<(std::ostream& out, const Field& field )
1078 cout << "Field " << field.getName() << " : " << endl ;
1079 out<< field.getField().retn()->getArray()->repr();