[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
|
From: |
Yves Renard |
|
Subject: |
[Getfem-commits] (no subject) |
|
Date: |
Sun, 6 May 2018 14:55:20 -0400 (EDT) |
branch: devel-yves-generic-assembly-modifs
commit 1258d6219fa49c2e3afed8583e6a072fa30d8062
Author: Yves Renard <address@hidden>
Date: Sun May 6 14:55:42 2018 +0200
Minor fixes
---
interface/tests/python/check_asm.py | 42 ++++++++++++++++++++++---
src/getfem_generic_assembly_compile_and_exec.cc | 15 ++++++---
src/getfem_generic_assembly_semantic.cc | 5 +--
3 files changed, 51 insertions(+), 11 deletions(-)
diff --git a/interface/tests/python/check_asm.py
b/interface/tests/python/check_asm.py
index ad8bf86..a3359d3 100644
--- a/interface/tests/python/check_asm.py
+++ b/interface/tests/python/check_asm.py
@@ -35,7 +35,7 @@ import os
NX = 4
m = gf.Mesh('triangles grid', np.arange(0,1+1./NX,1./NX),
np.arange(0,1+1./NX,1./NX)) # Structured mesh
-fem = gf.Fem('FEM_PK(2,1)')
+fem = gf.Fem('FEM_PK(2,1)')
mfu = gf.MeshFem(m, 1); mfu.set_fem(fem) # Lagrange P1 scalar fem
mim = gf.MeshIm(m, gf.Integ('IM_TRIANGLE(4)'))
mfv = gf.MeshFem(m, 3); mfv.set_fem(fem) # Lagrange P1 vector fem
@@ -52,10 +52,44 @@ md.set_variable('u', U)
md.add_fem_variable('v', mfv)
md.set_variable('v', V)
+# Simple test on the integral of u
+result = gf.asm('generic', mim, 0, "u", -1, md)
+if (abs(result-1) > 1e-8) : print "Bad value"; exit(1)
-print gf.asm('generic', mim, 0, "Def P(a):=a*(a'); Print(P(Grad_v))", -1, md)
-print gf.asm('generic', mim, 1, "Print(Grad_Test_u)", -1, md)
-print gf.asm('generic', mim, 0, "Def P(a):=a*(a'); Contract(P(Grad_v), 1, 2)",
-1, md)
+# Single contraction and comparison with Trace
+result1 = gf.asm('generic', mim, 0,
+ "Def P(a):=a*(a'); Contract(P(Grad_v), 1, 2)", -1, md)
+result2 = gf.asm('generic', mim, 0,
+ "Def P(a):=a*(a'); Trace(P(Grad_v))", -1, md)
+if (abs(result1-result2) > 1e-8) : print "Bad value"; exit(1)
+
+# Constant order 3 tensor contraction test
+result1 = gf.asm('generic', mim, 0,
+ "Contract([[[1,1],[2,2]],[[1,1],[2,2]]], 1, 2)", -1, md)
+result2 = np.array([3., 3.]);
+if (np.linalg.norm(result1-result2) > 1e-8) : print "Bad value"; exit(1)
+
+# Single contraction, comparison with "*"
+result1 = gf.asm('generic', mim, 0, "Contract(Grad_v, 2, Grad_u, 1)", -1, md)
+result2 = gf.asm('generic', mim, 0, "Grad_v * Grad_u", -1, md)
+if (np.linalg.norm(result1-result2) > 1e-8) : print "Bad value"; exit(1)
+
+# Double contraction order one expression, comparison with ":"
+result1 = gf.asm('generic', mim, 1,
+ "Contract(Grad_v, 1, 2, Grad_Test_v, 1, 2)", -1, md)
+result2 = gf.asm('generic', mim, 1, "Grad_v : Grad_Test_v", -1, md)
+if (np.linalg.norm(result1-result2) > 1e-8) : print "Bad value"; exit(1)
+
+# Double contraction order two expression, comparison with ":"
+result1 = gf.asm('generic', mim, 2,
+ "Contract(Grad_Test2_v, 1, 2, Grad_Test_v, 1, 2)", -1, md)
+result2 = gf.asm('generic', mim, 2, "Grad_Test2_v : Grad_Test_v", -1, md)
+if (np.linalg.norm(result1.full()-result2.full()) > 1e-8) :
+ print "Bad value"; exit(1)
+result1 = gf.asm('generic', mim, 2,
+ "Contract(Grad_Test_v, 2, 1, Grad_Test2_v, 2, 1)", -1, md)
+if (np.linalg.norm(result1.full()-result2.full()) > 1e-8) :
+ print "Bad value"; exit(1)
diff --git a/src/getfem_generic_assembly_compile_and_exec.cc
b/src/getfem_generic_assembly_compile_and_exec.cc
index a341d6d..af7d52a 100644
--- a/src/getfem_generic_assembly_compile_and_exec.cc
+++ b/src/getfem_generic_assembly_compile_and_exec.cc
@@ -1969,7 +1969,7 @@ namespace getfem {
virtual int exec() {
GA_DEBUG_INFO("Instruction: single contraction on a single tensor");
- size_type ii1 = tc1.size() / (nn*ii2*ii3);
+ size_type ii1 = tc1.size() / (nn*nn*ii2*ii3);
base_tensor::iterator it = t.begin();
for (size_type i = 0; i < ii3; ++i)
@@ -5923,7 +5923,7 @@ namespace getfem {
size_type kk = 0, ll = 1;
for (size_type i = 1; i < pnode->children.size(); ++i) {
if (i == ind[kk]) {
- ind[kk] = size_type(round(pnode->children[i]->tensor()[0]));
+ ind[kk] = size_type(round(pnode->children[i]->tensor()[0])-1);
indsize[kk] = pnode->children[ll]->tensor_proper_size(ind[kk]);
++kk;
} else ll = i;
@@ -5963,10 +5963,11 @@ namespace getfem {
}
else if (pnode->children.size() == 7) {
// Particular cases should be detected (ii2=ii3=1 in particular).
- GMM_ASSERT1(false, "to be done !");
size_type i1 = ind[0], i2 = ind[1], i3 = ind[2], i4 = ind[3];
size_type nn1 = indsize[0], nn2 = indsize[1];
size_type ii1 = 1, ii2 = 1, ii3 = 1, ii4 = 1, ii5 = 1, ii6 = 1;
+ if (i1 > i2)
+ { std::swap(i1, i2); std::swap(i3, i4); std::swap(nn1, nn2); }
for (size_type i = 0; i < child1->tensor_order(); ++i) {
if (i < i1) ii1 *= child1->tensor_proper_size(i);
if (i > i1 && i < i2) ii2 *= child1->tensor_proper_size(i);
@@ -5978,8 +5979,6 @@ namespace getfem {
ii5 *= child2->tensor_proper_size(i);
if (i > i3 && i > i4) ii6 *= child2->tensor_proper_size(i);
}
- if (i1 > i2)
- { std::swap(i1, i2); std::swap(i3, i4); std::swap(nn1, nn2); }
if (child1->test_function_type==1 && child2->test_function_type==2)
pgai = std::make_shared<ga_instruction_contract_2_2_rev>
(pnode->tensor(), child1->tensor(), child2->tensor(),
@@ -6331,6 +6330,9 @@ namespace getfem {
break;
case 1:
{
+ GMM_ASSERT1(root->tensor_proper_size() == 1,
+ "Invalid vector or tensor quantity. An order 1 "
+ "weak form has to be a scalar quantity");
const mesh_fem *mf=workspace.associated_mf(root->name_test1);
const mesh_fem **mfg = 0;
add_interval_to_gis(workspace, root->name_test1, gis);
@@ -6370,6 +6372,9 @@ namespace getfem {
break;
case 2:
{
+ GMM_ASSERT1(root->tensor_proper_size() == 1,
+ "Invalid vector or tensor quantity. An order 2 "
+ "weak form has to be a scalar quantity");
const mesh_fem
*mf1=workspace.associated_mf(root->name_test1);
const mesh_fem
*mf2=workspace.associated_mf(root->name_test2);
const mesh_fem **mfg1 = 0, **mfg2 = 0;
diff --git a/src/getfem_generic_assembly_semantic.cc
b/src/getfem_generic_assembly_semantic.cc
index 6bd1c8d..07340cb 100644
--- a/src/getfem_generic_assembly_semantic.cc
+++ b/src/getfem_generic_assembly_semantic.cc
@@ -1800,12 +1800,13 @@ namespace getfem {
"Invalid parameter for Contract. "
"Should be an index number.");
ind[kk] = size_type(round(pnode->children[i]->tensor()[0]));
- indsize[kk] = pnode->children[ll]->tensor_proper_size(ind[kk]);
order = pnode->children[ll]->tensor_order();
if (ind[kk] < 1 || ind[kk] > order)
ga_throw_error(pnode->children[i]->expr, pnode->children[i]->pos,
"Parameter out of range for Contract (should be "
"between 1 and " << order << ")");
+ ind[kk]--;
+ indsize[kk] = pnode->children[ll]->tensor_proper_size(ind[kk]);
if (kk >= ind.size()/2 && indsize[kk] != indsize[kk-ind.size()/2])
ga_throw_error(child0->expr, child0->pos,
"Invalid parameters for Contract. Cannot "
@@ -1827,7 +1828,7 @@ namespace getfem {
size_type i1 = ind[0], i2 = ind[1];
if (i1 == i2)
ga_throw_error(child0->expr, child0->pos,
- "Invalid parameters for Contract. Repeated index");
+ "Invalid parameters for Contract. Repeated index.");
mi.resize(0);
for (size_type i = 0; i < pnode->nb_test_functions(); ++i)