From 556c94b24e137571a459169d28552ee24b82ddfe Mon Sep 17 00:00:00 2001
From: Harry Cousin <harry.cousin@inrae.fr>
Date: Thu, 3 Oct 2024 15:30:54 +0200
Subject: [PATCH 1/8] add new argument int* jMap in addMatMapij
---
trunk/mseMatTools/src/mseMatTools.c | 8 ++++----
trunk/mseMatTools/src/mseMatTools.h | 2 +-
2 files changed, 5 insertions(+), 5 deletions(-)
diff --git a/trunk/mseMatTools/src/mseMatTools.c b/trunk/mseMatTools/src/mseMatTools.c
index edd51c6..f3e952e 100644
--- a/trunk/mseMatTools/src/mseMatTools.c
+++ b/trunk/mseMatTools/src/mseMatTools.c
@@ -56,7 +56,7 @@ void assembleMat(Mat matOut, Mat matIn,int NOffset, int MOffset){
/*
*#pMatOut = pMatOut + pMatIn[iMap[i],iMap[j]]
**/
-void addMatMapij(Mat matOut, double alpha, Mat matIn, int *iMap, int iLength, int iMask, int jMask){
+void addMatMapij(Mat matOut, double alpha, Mat matIn, int *iMap, int *jMap, int iLength, int iMask, int jMask){
PetscInt NO,MO;
MatGetSize(matOut,&NO,&MO);
PetscInt NI,MI,i,ii,iOffseted;
@@ -71,11 +71,11 @@ void addMatMapij(Mat matOut, double alpha, Mat matIn, int *iMap, int iLength, in
#ifdef mseMatToolsS_VERBOSE
printf("mseMatToolsS addMatMapij No=%i Mo=%i Ni=%i Mi=%i.\n",NO,MO,NI,MI);
if (iMask && jMask)
- printf(" doing matOut[iMap[i],iMap[j]]+=matIn[i,j]\n");
+ printf(" doing matOut[iMap[i],jMap[j]]+=matIn[i,j]\n");
else if (iMask)
printf(" doing matOut[iMap[i],j]+=matIn[i,j]\n");
else if (jMask)
- printf(" doing matOut[i,iMap[j]]+=matIn[i,j]\n");
+ printf(" doing matOut[i,jMap[j]]+=matIn[i,j]\n");
#endif
for (i = 0; i < NI; i++) {
@@ -89,7 +89,7 @@ void addMatMapij(Mat matOut, double alpha, Mat matIn, int *iMap, int iLength, in
printf("mseMatToolsS addMatMapij row %i, with %i nc\n",i,nc);
#endif
for (ii=0;ii<nc;ii++){
- ajOffseted[ii]=jMask?iMap[aj[ii]]:aj[ii];
+ ajOffseted[ii]=jMask?jMap[aj[ii]]:aj[ii];
copyaa[ii]=alpha* copyaa[ii];
}
iOffseted=iMask?iMap[i]:i;
diff --git a/trunk/mseMatTools/src/mseMatTools.h b/trunk/mseMatTools/src/mseMatTools.h
index c50e4ad..6d4f5ac 100644
--- a/trunk/mseMatTools/src/mseMatTools.h
+++ b/trunk/mseMatTools/src/mseMatTools.h
@@ -11,7 +11,7 @@ void assembleVec(Vec pVecOut,Vec pVecIn,int Offset);
*#pMatOut = pMatOut + pMatIn[iMap[i],iMap[j]]
* memo: numpy.i export array in (int *iMap, int iLength)
**/
-void addMatMapij(Mat pMatOut, double alpha, Mat pMatIn, int *iMap, int iLength, int iMask, int jMask);
+void addMatMapij(Mat pMatOut, double alpha, Mat pMatIn, int *iMap, int* jMap, int iLength, int iMask, int jMask);
void printMat(Mat aMat);
void saveVec(char* fileName,double *dVec, int dLength);
void loadVec(char* fileName,double *dVec, int dLength);
--
GitLab
From 1a51d50a7de77efccaabfc675d82321d48771fde Mon Sep 17 00:00:00 2001
From: Harry Cousin <harry.cousin@inrae.fr>
Date: Fri, 4 Oct 2024 11:31:05 +0200
Subject: [PATCH 2/8] TO FIX: modif to add new argument, study3 doesn't work
---
trunk/mseMatTools/src/mseMatTools.c | 2 +-
trunk/mseMatTools/src/mseMatTools.h | 2 +-
trunk/src/mse/INTERACTIONS/interaction.py | 6 +++---
trunk/src/mse/linAlgTools.py | 4 ++--
4 files changed, 7 insertions(+), 7 deletions(-)
diff --git a/trunk/mseMatTools/src/mseMatTools.c b/trunk/mseMatTools/src/mseMatTools.c
index f3e952e..64c91d2 100644
--- a/trunk/mseMatTools/src/mseMatTools.c
+++ b/trunk/mseMatTools/src/mseMatTools.c
@@ -56,7 +56,7 @@ void assembleMat(Mat matOut, Mat matIn,int NOffset, int MOffset){
/*
*#pMatOut = pMatOut + pMatIn[iMap[i],iMap[j]]
**/
-void addMatMapij(Mat matOut, double alpha, Mat matIn, int *iMap, int *jMap, int iLength, int iMask, int jMask){
+void addMatMapij(Mat matOut, double alpha, Mat matIn, int *iMap, int iLength, int *jMap, int jLength, int iMask, int jMask){
PetscInt NO,MO;
MatGetSize(matOut,&NO,&MO);
PetscInt NI,MI,i,ii,iOffseted;
diff --git a/trunk/mseMatTools/src/mseMatTools.h b/trunk/mseMatTools/src/mseMatTools.h
index 6d4f5ac..b8fd108 100644
--- a/trunk/mseMatTools/src/mseMatTools.h
+++ b/trunk/mseMatTools/src/mseMatTools.h
@@ -11,7 +11,7 @@ void assembleVec(Vec pVecOut,Vec pVecIn,int Offset);
*#pMatOut = pMatOut + pMatIn[iMap[i],iMap[j]]
* memo: numpy.i export array in (int *iMap, int iLength)
**/
-void addMatMapij(Mat pMatOut, double alpha, Mat pMatIn, int *iMap, int* jMap, int iLength, int iMask, int jMask);
+void addMatMapij(Mat pMatOut, double alpha, Mat pMatIn, int *iMap, int iLength, int* jMap, int jLength, int iMask, int jMask);
void printMat(Mat aMat);
void saveVec(char* fileName,double *dVec, int dLength);
void loadVec(char* fileName,double *dVec, int dLength);
diff --git a/trunk/src/mse/INTERACTIONS/interaction.py b/trunk/src/mse/INTERACTIONS/interaction.py
index 1d76fed..658ee35 100644
--- a/trunk/src/mse/INTERACTIONS/interaction.py
+++ b/trunk/src/mse/INTERACTIONS/interaction.py
@@ -127,10 +127,10 @@ class interaction2D1D(interaction):
self.addM2D1D.setDim(self.dS2D.nbDofsPerComp,self.dS1D.nbDofsPerComp)
self.addM1D2D=matrix()
self.addM1D2D.setDim(self.dS1D.nbDofsPerComp,self.dS2D.nbDofsPerComp)
- self.addM2D.addMatMapij(1.0,self.dS1D.massMat,self.index1DTo2D,1,1)
+ self.addM2D.addMatMapij(1.0,self.dS1D.massMat,self.index1DTo2D,self.index1DTo2D,1,1)
self.addM1D.axpy(1.0,self.dS1D.massMat)
- self.addM2D1D.addMatMapij(-1.0,self.dS1D.massMat,self.index1DTo2D,1,0)
- self.addM1D2D.addMatMapij(-1.0,self.dS1D.massMat,self.index1DTo2D,0,1)
+ self.addM2D1D.addMatMapij(-1.0,self.dS1D.massMat,self.index1DTo2D,self.index1DTo2D,1,0)
+ self.addM1D2D.addMatMapij(-1.0,self.dS1D.massMat,self.index1DTo2D,self.index1DTo2D,0,1)
for i in range(self.dS1.nbComps):
if (self.compsInvolved[i]):
#self.addM1=self.addM2D
diff --git a/trunk/src/mse/linAlgTools.py b/trunk/src/mse/linAlgTools.py
index ddd6f03..50da82b 100644
--- a/trunk/src/mse/linAlgTools.py
+++ b/trunk/src/mse/linAlgTools.py
@@ -74,8 +74,8 @@ class matrix:
# if Jmask and not Imask:
# self[i,iMap[j]] = self[i,iMap[j]] + alpha*aMat[i,j] ou i,j index de aMat
- def addMatMapij(self,alpha,aMat,iMap,Imask,Jmask):
- MT.addMatMapij(self.petscMat,alpha,aMat.petscMat,iMap,Imask,Jmask)
+ def addMatMapij(self,alpha,aMat,iMap,jMap,Imask,Jmask):
+ MT.addMatMapij(self.petscMat,alpha,aMat.petscMat,iMap,jMap,Imask,Jmask)
class vector:
--
GitLab
From 3cf5e22b21f61684601fbb1fab1abd5899e0f91d Mon Sep 17 00:00:00 2001
From: Harry Cousin <harry.cousin@inrae.fr>
Date: Fri, 4 Oct 2024 17:58:15 +0200
Subject: [PATCH 3/8] modif in file mseMatTools.i.in
---
trunk/mseMatTools/Frontend/mseMatTools.i.in | 1 +
1 file changed, 1 insertion(+)
diff --git a/trunk/mseMatTools/Frontend/mseMatTools.i.in b/trunk/mseMatTools/Frontend/mseMatTools.i.in
index c4c9ccc..dec4960 100644
--- a/trunk/mseMatTools/Frontend/mseMatTools.i.in
+++ b/trunk/mseMatTools/Frontend/mseMatTools.i.in
@@ -11,6 +11,7 @@
%}
%apply (int *IN_ARRAY1, int DIM1) {(int* iMap, int iLength)}
+%apply (int *IN_ARRAY1, int DIM1) {(int* jMap, int jLength)}
%apply (double *IN_ARRAY1, int DIM1) {(double *dVec, int dLength)}
%apply (double *IN_ARRAY1, int DIM1) {(double *dVec1, int dLength1)}
%apply (double *IN_ARRAY1, int DIM1) {(double *dVec2, int dLength2)}
--
GitLab
From 29ba617d9190799849ce093433d9c75810689b99 Mon Sep 17 00:00:00 2001
From: Harry Cousin <harry.cousin@inrae.fr>
Date: Tue, 8 Oct 2024 09:30:08 +0200
Subject: [PATCH 4/8] add class interaction1D1D, without test
---
trunk/src/mse/INTERACTIONS/interaction.py | 31 +++++++++++++++++++++++
1 file changed, 31 insertions(+)
diff --git a/trunk/src/mse/INTERACTIONS/interaction.py b/trunk/src/mse/INTERACTIONS/interaction.py
index 658ee35..f793076 100644
--- a/trunk/src/mse/INTERACTIONS/interaction.py
+++ b/trunk/src/mse/INTERACTIONS/interaction.py
@@ -141,3 +141,34 @@ class interaction2D1D(interaction):
self.addM12.assembleMat(self.addM2D1D,i*self.dS2.nbDofsPerComp,i*self.dS1.nbDofsPerComp)
#self.addM21=self.addM1D2D
self.addM21.assembleMat(self.addM1D2D,i*self.dS1.nbDofsPerComp,i*self.dS2.nbDofsPerComp)
+###
+# CLASS INTERACTION 1D1D
+###
+class interaction1D1D(interaction):
+ def __init__(self, aStudy, aDS1d1, aDS1d2, aCompsInvolved=None) -> None:
+ super().init__init__(aStudy,aDS1d1,aDS1d2,aCompsInvolved)
+ self.dS1D1=aDS1d1
+ self.dS1D2=aDS1d2
+ if (self.dS1D1.geoDim != 1 or self.dS1D1.geoDim != 1):
+ print("ERROR in interaction2D1D, wrong dimension " +str(self.dS2D.aGeoDim)+" "+self.dS1D.aGeoDim)
+ exit()
+ # ?
+ self.c1d2d=1
+ self.c2d1D=1
+ #
+ self.indexG1ToG2=np.zeros((self.dS1D1.nbDofsPerComp,),dtype=np.int32)
+ self.indexG2ToG1=np.zeros((self.dS1D2.nbDofsPerComp,),dtype=np.int32)# dS1D1.nbDof == dS1D2.nbDof
+ buildMapG1ToG2(self.dS1D1,self.dS1D2,self.indexG1ToG2)
+ buildMapG2ToG1(self.dS1D2,self.dS1D1,self.indexG2ToG1)
+ self.addM2=matrix()
+ self.addM2.setDim(self.dS1D2.nbDofsPerComp,self.dS1D2.nbDofsPerComp)
+ self.addM1=matrix()
+ self.addM1.setDim(self.dS1D1.nbDofsPerComp,self.dS1D1.nbDofsPerComp)
+ self.addM21=matrix()
+ self.addM21.setDim(self.dS1D2.nbDofsPerComp,self.dS1D1.nbDofsPerComp)
+ self.addM12=matrix()
+ self.addM12.setDim(self.dS1D1.nbDofsPerComp,self.dS1D1.nbDofsPerComp)
+ self.addM2.axpy(1.0,self.dS1D2.massMat)
+ self.addM1.axpy(1.0,self.dS1D1.massMat)
+ self.addM21.addMatMapij(-1.0,self.dS1D.massMat,self.index2DTo1D,self.index2DTo1D,0,1)
+ self.addM12.addMatMapij(-1.0,self.dS1D.massMat,self.index1DTo2D,self.index1DTo2D,0,1)
--
GitLab
From 2f540686a6b05b76c59a14186f8bdf09d4fa1efd Mon Sep 17 00:00:00 2001
From: Harry Cousin <harry.cousin@inrae.fr>
Date: Tue, 8 Oct 2024 16:55:04 +0200
Subject: [PATCH 5/8] patch interaction 1D1D
---
trunk/src/mse/INTERACTIONS/interaction.py | 10 +++++-----
1 file changed, 5 insertions(+), 5 deletions(-)
diff --git a/trunk/src/mse/INTERACTIONS/interaction.py b/trunk/src/mse/INTERACTIONS/interaction.py
index f793076..3639808 100644
--- a/trunk/src/mse/INTERACTIONS/interaction.py
+++ b/trunk/src/mse/INTERACTIONS/interaction.py
@@ -146,7 +146,7 @@ class interaction2D1D(interaction):
###
class interaction1D1D(interaction):
def __init__(self, aStudy, aDS1d1, aDS1d2, aCompsInvolved=None) -> None:
- super().init__init__(aStudy,aDS1d1,aDS1d2,aCompsInvolved)
+ super().__init__(aStudy,aDS1d1,aDS1d2,aCompsInvolved)
self.dS1D1=aDS1d1
self.dS1D2=aDS1d2
if (self.dS1D1.geoDim != 1 or self.dS1D1.geoDim != 1):
@@ -158,8 +158,8 @@ class interaction1D1D(interaction):
#
self.indexG1ToG2=np.zeros((self.dS1D1.nbDofsPerComp,),dtype=np.int32)
self.indexG2ToG1=np.zeros((self.dS1D2.nbDofsPerComp,),dtype=np.int32)# dS1D1.nbDof == dS1D2.nbDof
- buildMapG1ToG2(self.dS1D1,self.dS1D2,self.indexG1ToG2)
- buildMapG2ToG1(self.dS1D2,self.dS1D1,self.indexG2ToG1)
+ buildMap1DTo2D(self.dS1D1,self.dS1D2,self.indexG1ToG2)
+ buildMap1DTo2D(self.dS1D2,self.dS1D1,self.indexG2ToG1)
self.addM2=matrix()
self.addM2.setDim(self.dS1D2.nbDofsPerComp,self.dS1D2.nbDofsPerComp)
self.addM1=matrix()
@@ -170,5 +170,5 @@ class interaction1D1D(interaction):
self.addM12.setDim(self.dS1D1.nbDofsPerComp,self.dS1D1.nbDofsPerComp)
self.addM2.axpy(1.0,self.dS1D2.massMat)
self.addM1.axpy(1.0,self.dS1D1.massMat)
- self.addM21.addMatMapij(-1.0,self.dS1D.massMat,self.index2DTo1D,self.index2DTo1D,0,1)
- self.addM12.addMatMapij(-1.0,self.dS1D.massMat,self.index1DTo2D,self.index1DTo2D,0,1)
+ self.addM21.addMatMapij(-1.0,self.dS1D2.massMat,self.indexG2ToG1,self.indexG2ToG1,0,1)
+ self.addM12.addMatMapij(-1.0,self.dS1D1.massMat,self.indexG1ToG2,self.indexG1ToG2,0,1)
--
GitLab
From 59e0aea89619416b1d477ec4bae61f1560146a83 Mon Sep 17 00:00:00 2001
From: Harry Cousin <harry.cousin@inrae.fr>
Date: Tue, 8 Oct 2024 16:56:31 +0200
Subject: [PATCH 6/8] new study to test class interaction1D1D + new file
sourceModel affineGrowth with param
---
trunk/BD_MODELS/affineGrowthParam.txt | 1 +
trunk/DEMOS/study14.py | 106 ++++++++++++++++++++++++++
2 files changed, 107 insertions(+)
create mode 100644 trunk/BD_MODELS/affineGrowthParam.txt
create mode 100644 trunk/DEMOS/study14.py
diff --git a/trunk/BD_MODELS/affineGrowthParam.txt b/trunk/BD_MODELS/affineGrowthParam.txt
new file mode 100644
index 0000000..7c268fa
--- /dev/null
+++ b/trunk/BD_MODELS/affineGrowthParam.txt
@@ -0,0 +1 @@
+Fa=p1-a
diff --git a/trunk/DEMOS/study14.py b/trunk/DEMOS/study14.py
new file mode 100644
index 0000000..bae64fb
--- /dev/null
+++ b/trunk/DEMOS/study14.py
@@ -0,0 +1,106 @@
+from mse.study import study
+from mse.dynamicalSystem import dynamicalSystem,computeMass
+from mse.DISPERSIONS.diffusionModel import diffusionModel
+from mse.MODELS.models import txtModel
+from mse.system import system
+from mse.simulator import simulator
+from mse.timeDoOneStep import timeDoOneStep
+from mse.TIMESCHEMES.timeScheme import implicitLinearTimeScheme,implicitNonLinearTimeScheme
+from mse.INTERACTIONS.interaction import interaction2D1D, interaction1D1D
+from mse.toolsEnv import computerEnv,MSEPrint
+from mse.explorTrajectory import explorTrajectory
+from os import mkdir
+from os.path import isdir
+
+
+nbComps=1
+WITH_AFFINE_GROWTH=1
+testNLTS=0
+studName="STUDY14_DIFFUSION"
+if (WITH_AFFINE_GROWTH):
+ studName=studName+"_AND_AFFINE_GROWTH"
+nbParam=1
+#a study
+aStudy=study(studName, nbParam)
+# the study zone
+aStudy.setGeometry("SIMPLE2",0.4)
+aStudy.setParam([0.5])
+
+# a system
+aSystem=system(aStudy)
+# dynamical system
+# First DS
+ds2d1=dynamicalSystem(aStudy,2,"surface1.xdmf",nbComps=nbComps)
+# Second DS
+ds2d2=dynamicalSystem(aStudy,2,"surface2.xdmf",nbComps=nbComps)
+# edge
+ds1d1=dynamicalSystem(aStudy,1,"edge2.xdmf",nbComps=nbComps)
+ds1d2=dynamicalSystem(aStudy,1,"edge2.xdmf",nbComps=nbComps)
+y0=-0.45
+coefInit = 1
+ds2d1.initialState[0]="np.exp(-"+str(coefInit)+"*(x[0]*0.0+(x[1]-"+str(y0)+")**2))"
+ds2d2.initialState[0]="0.0*"+ds2d2.initialState[0]
+ds1d1.initialState[0]="np.exp(-"+str(coefInit)+"*(x[0]*0.0+(x[1]-"+str(y0)+")**2))"
+ds1d2.initialState[0]="0.0*"+ds2d2.initialState[0]
+if (nbComps>1):
+ ds2d1.initialState[1]="0+0.0*"+ds2d1.initialState[0]
+ ds2d2.initialState[1]="0.0+0.0*"+ds2d1.initialState[0]
+ ds1d1.initialState[1]="0.0+0.0*"+ds2d1.initialState[0]
+ ds1d2.initialState[1]="0.0+0.0*"+ds2d1.initialState[0]
+#
+
+#add diffusion motion
+diffCoefOmega = 0.2
+diffusionModel(ds2d1).setDiffCoef(diffCoefOmega)
+diffusionModel(ds2d2).setDiffCoef(1.0)
+diffusionModel(ds1d1).setDiffCoef(3.)
+diffusionModel(ds1d2).setDiffCoef(1.)
+if (WITH_AFFINE_GROWTH):
+ m=txtModel(ds2d1,"affineGrowthParam")
+ m.computeMatOnce=1
+ m.computeRhsOnce=1
+ m1d1=txtModel(ds1d1,"affineGrowthParam")
+ m1d1.computeMatOnce=1
+ m1d1.computeRhsOnce=1
+ m1d2=txtModel(ds1d2,"affineGrowthParam")
+ m1d2.computeMatOnce=1
+ m1d2.computeRhsOnce=1
+
+aInter2D1D1=interaction2D1D(aStudy,ds2d1,ds1d1)
+aInter2D1D2=interaction2D1D(aStudy,ds2d2,ds1d2)
+
+aInter1D1D=interaction1D1D(aStudy, ds1d1,ds1d2)
+
+
+
+
+aTS=None
+# select the time algorithm
+if (testNLTS):
+ aTS=implicitNonLinearTimeScheme()
+else:
+ aTS=implicitLinearTimeScheme()
+#simu in t \in [0,1]
+simulator(0,1,aStudy)
+#based on fix time step 0.05
+aTStep=timeDoOneStep(aStudy,aTS,0.05)
+#run the simulation
+aStudy.simulator.doSimu()
+
+
+#for visu, export simulation steps to vtk
+aStudy.simulator.exportSimuToVtk()
+
+#loop on simu to compute mass
+import numpy as np
+for ds in aStudy.system.dynamicalSystems:
+ mass=np.zeros((ds.nbComps),dtype=np.double)
+ aExplorTraj=explorTrajectory(aStudy,aStudy.simuDir)
+ aExplorTraj.rewind()
+ while(aExplorTraj.replay()):
+ computeMass(ds,mass)
+ print("mass "+str(aExplorTraj.curStep)+" = ",end="")
+ for m in mass:
+ print(m,end=", ")
+ print("")
+aStudy.end()
--
GitLab
From c67dd6ed2bc851191d29bbc46f9bb4c671f812a6 Mon Sep 17 00:00:00 2001
From: Harry Cousin <harry.cousin@inrae.fr>
Date: Mon, 14 Oct 2024 10:45:49 +0200
Subject: [PATCH 7/8] option in study to add or not interaciont 1d1d + begin
creation of class interaction2D2D
---
trunk/DEMOS/study14.py | 6 +++-
trunk/src/mse/INTERACTIONS/interaction.py | 34 +++++++++++++++++++++++
2 files changed, 39 insertions(+), 1 deletion(-)
diff --git a/trunk/DEMOS/study14.py b/trunk/DEMOS/study14.py
index bae64fb..dafd951 100644
--- a/trunk/DEMOS/study14.py
+++ b/trunk/DEMOS/study14.py
@@ -15,10 +15,13 @@ from os.path import isdir
nbComps=1
WITH_AFFINE_GROWTH=1
+WITH_INT_1D1D=1
testNLTS=0
studName="STUDY14_DIFFUSION"
if (WITH_AFFINE_GROWTH):
studName=studName+"_AND_AFFINE_GROWTH"
+if (WITH_INT_1D1D):
+ studName=studName+"_AND_INT_1D1D"
nbParam=1
#a study
aStudy=study(studName, nbParam)
@@ -69,7 +72,8 @@ if (WITH_AFFINE_GROWTH):
aInter2D1D1=interaction2D1D(aStudy,ds2d1,ds1d1)
aInter2D1D2=interaction2D1D(aStudy,ds2d2,ds1d2)
-aInter1D1D=interaction1D1D(aStudy, ds1d1,ds1d2)
+if (WITH_AFFINE_GROWTH):
+ aInter1D1D=interaction1D1D(aStudy, ds1d1,ds1d2)
diff --git a/trunk/src/mse/INTERACTIONS/interaction.py b/trunk/src/mse/INTERACTIONS/interaction.py
index 3639808..ed5315b 100644
--- a/trunk/src/mse/INTERACTIONS/interaction.py
+++ b/trunk/src/mse/INTERACTIONS/interaction.py
@@ -172,3 +172,37 @@ class interaction1D1D(interaction):
self.addM1.axpy(1.0,self.dS1D1.massMat)
self.addM21.addMatMapij(-1.0,self.dS1D2.massMat,self.indexG2ToG1,self.indexG2ToG1,0,1)
self.addM12.addMatMapij(-1.0,self.dS1D1.massMat,self.indexG1ToG2,self.indexG1ToG2,0,1)
+###
+###
+# CLASS INTERACTION 2D2D
+###
+#class interaction2D2D(interaction):
+# def __init__(self, aStudy, aDS2d1, aDS2d2, aCompsInvolved=None) -> None:
+# super().__init__(aStudy,aDS1d1,aDS1d2,aCompsInvolved)
+# self.dS2D1=aDS2d1
+# self.dS2D2=aDS2d2
+# if (self.dS1D1.geoDim != 2 or self.dS1D1.geoDim != 2):
+# print("ERROR in interaction2D1D, wrong dimension " +str(self.dS2D.aGeoDim)+" "+self.dS1D.aGeoDim)
+# exit()
+# # ?
+# self.c1d2d=1
+# self.c2d1D=1
+# #
+# #TODO: change all
+# self.indexG1ToG2=np.zeros((self.dS1D1.nbDofsPerComp,),dtype=np.int32)
+# self.indexG2ToG1=np.zeros((self.dS1D2.nbDofsPerComp,),dtype=np.int32)# dS1D1.nbDof == dS1D2.nbDof
+# buildMap1DTo2D(self.dS1D1,self.dS1D2,self.indexG1ToG2)
+# buildMap1DTo2D(self.dS1D2,self.dS1D1,self.indexG2ToG1)
+# self.addM2=matrix()
+# self.addM2.setDim(self.dS1D2.nbDofsPerComp,self.dS1D2.nbDofsPerComp)
+# self.addM1=matrix()
+# self.addM1.setDim(self.dS1D1.nbDofsPerComp,self.dS1D1.nbDofsPerComp)
+# self.addM21=matrix()
+# self.addM21.setDim(self.dS1D2.nbDofsPerComp,self.dS1D1.nbDofsPerComp)
+# self.addM12=matrix()
+# self.addM12.setDim(self.dS1D1.nbDofsPerComp,self.dS1D1.nbDofsPerComp)
+# self.addM2.axpy(1.0,self.dS1D2.massMat)
+# self.addM1.axpy(1.0,self.dS1D1.massMat)
+# self.addM21.addMatMapij(-1.0,self.dS1D2.massMat,self.indexG2ToG1,self.indexG2ToG1,0,1)
+# self.addM12.addMatMapij(-1.0,self.dS1D1.massMat,self.indexG1ToG2,self.indexG1ToG2,0,1)
+# # END TODO
--
GitLab
From 16d817126f7e4a7dda2fd355e6ac4514d4cd92d2 Mon Sep 17 00:00:00 2001
From: Harry Cousin <harry.cousin@inrae.fr>
Date: Mon, 14 Oct 2024 14:53:31 +0200
Subject: [PATCH 8/8] patch study, add total compute mass
---
trunk/DEMOS/study14.py | 41 ++++++++++++++++++++++-------------------
1 file changed, 22 insertions(+), 19 deletions(-)
diff --git a/trunk/DEMOS/study14.py b/trunk/DEMOS/study14.py
index dafd951..87abd66 100644
--- a/trunk/DEMOS/study14.py
+++ b/trunk/DEMOS/study14.py
@@ -14,20 +14,19 @@ from os.path import isdir
nbComps=1
-WITH_AFFINE_GROWTH=1
+WITH_LOGISTIC=1
WITH_INT_1D1D=1
-testNLTS=0
+testNLTS=1
studName="STUDY14_DIFFUSION"
-if (WITH_AFFINE_GROWTH):
- studName=studName+"_AND_AFFINE_GROWTH"
+if (WITH_LOGISTIC):
+ studName=studName+"_AND_LOGISTIC"
if (WITH_INT_1D1D):
studName=studName+"_AND_INT_1D1D"
-nbParam=1
+nbParam=0
#a study
aStudy=study(studName, nbParam)
# the study zone
aStudy.setGeometry("SIMPLE2",0.4)
-aStudy.setParam([0.5])
# a system
aSystem=system(aStudy)
@@ -43,8 +42,8 @@ y0=-0.45
coefInit = 1
ds2d1.initialState[0]="np.exp(-"+str(coefInit)+"*(x[0]*0.0+(x[1]-"+str(y0)+")**2))"
ds2d2.initialState[0]="0.0*"+ds2d2.initialState[0]
-ds1d1.initialState[0]="np.exp(-"+str(coefInit)+"*(x[0]*0.0+(x[1]-"+str(y0)+")**2))"
-ds1d2.initialState[0]="0.0*"+ds2d2.initialState[0]
+ds1d1.initialState[0]="0.0*"+ds1d1.initialState[0]
+ds1d2.initialState[0]="0.0*"+ds1d2.initialState[0]
if (nbComps>1):
ds2d1.initialState[1]="0+0.0*"+ds2d1.initialState[0]
ds2d2.initialState[1]="0.0+0.0*"+ds2d1.initialState[0]
@@ -58,21 +57,13 @@ diffusionModel(ds2d1).setDiffCoef(diffCoefOmega)
diffusionModel(ds2d2).setDiffCoef(1.0)
diffusionModel(ds1d1).setDiffCoef(3.)
diffusionModel(ds1d2).setDiffCoef(1.)
-if (WITH_AFFINE_GROWTH):
- m=txtModel(ds2d1,"affineGrowthParam")
- m.computeMatOnce=1
- m.computeRhsOnce=1
- m1d1=txtModel(ds1d1,"affineGrowthParam")
- m1d1.computeMatOnce=1
- m1d1.computeRhsOnce=1
- m1d2=txtModel(ds1d2,"affineGrowthParam")
- m1d2.computeMatOnce=1
- m1d2.computeRhsOnce=1
+if (WITH_LOGISTIC):
+ m2=txtModel(ds2d1,"logistic")
aInter2D1D1=interaction2D1D(aStudy,ds2d1,ds1d1)
aInter2D1D2=interaction2D1D(aStudy,ds2d2,ds1d2)
-if (WITH_AFFINE_GROWTH):
+if (WITH_INT_1D1D):
aInter1D1D=interaction1D1D(aStudy, ds1d1,ds1d2)
@@ -107,4 +98,16 @@ for ds in aStudy.system.dynamicalSystems:
for m in mass:
print(m,end=", ")
print("")
+
+mass=np.zeros((ds.nbComps),dtype=np.double)
+aExplorTraj=explorTrajectory(aStudy,aStudy.simuDir)
+aExplorTraj.rewind()
+while(aExplorTraj.replay()):
+ TotalMass=np.zeros((ds.nbComps),dtype=np.double)
+ for ds in aStudy.system.dynamicalSystems:
+ computeMass(ds,mass)
+ TotalMass=TotalMass+mass
+ print("total mass "+str(aExplorTraj.curStep)+" = ",end="")
+ print(TotalMass,end=", ")
+ print("")
aStudy.end()
--
GitLab